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SUMMARY 


Results of the present theory on the normal approach elastohydro- 
dynamic problem show that: 

1. The features of pressure and deformation profiles during the 
early stages of the normal approach agree well with those 
obtained in Ref. 4, which neglects the influence of the 
local approach velocity. The steepness of the pressure 
gradient at the center is strongly dependent upon the 
product of the pressure-viscosity coefficient and the center 
pressure. This strong dependence is removed if a smaller 
pressure-viscosity coefficient is used at high pressures. 

2. During final stages of the normal approach, present theory 
yields considerably different results from those in Ref. 4. 
The local approach velocity at the edge of the contact 
region becomes far greater than the center approach velocity, 
and finally entraps a pocket of the lubricant at the center 
of the contact. Both the deformation and pressure profiles 
never converge to the dry contact Hertzian distribution. 

3. For a normal approach process under a constant load, the max- 

imum center pressure can exceed that of the maximum Hertzian 
pressure depending upon the pressure-viscosity coefficient. 

By introducing the composite-exponential model for the 
pressure-viscosity dependence, the maximum center pressure 
is much reduced. 

4. The inclusion of the lubricant compressibility in the analysis 
gives arise to a slightly higher load than the incompressible 
solution. 

iv 
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CHAPTER 1 - INTRODUCTION 

1.1 Introduction 

Whenever any two lubricated contacts approach each other along 
their common normal under a heavy load, highly localized pressures are 
generated by the squeeze film action within the conjunction. The deter- 
mination of the pressure distribution due to the squeeze action consider- 
ing the surface deformation is known as the normal approach problem in 
el as tohydrodynamic (EHD) lubrication. 

The squeeze-film action occurs frequently in many machine components 
such as gear teeth contacts, cams, and rolling element bearings during 
transient loadings. The normal approach problem has a special signifi- 
cance in the so-called partial EHD contacts in which the asperity heights 
approach the same order of magnitude as the film thickness. Under these 
conditions, the entering of any asperity into the conjunction zone is equi- 
valent to the squeeze-film EHD problem between a contacting body and a 
flat plate. 

Mathematically, the normal approach problem differs considerably 
from the conventional rolling and sliding EHD theories [l,2,3]. For 
the rolling problem, the pressure and film distributions are steady- 
state: whereas for squeeze-film problem they are time- dependent and 
must be obtained by solving the transient Reynolds equation coupled 
with the elasticity equation. Because the pressure gradient varies 
inversely with the third power of the film thickness and the viscosity 
for most lubricants varies exponentially with pressure, the two coupled 
equations are highly nonlinear. So far, no analytical solution has 
been found for these equations. 
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In 1961, Christensen [4] introduced the first numerical solution 
to the present EHD problem for an incompressible lubricant with an 
exponentially varying viscosity. In his solution, he has neglected 
the squeeze-film action due to the change of deformation. This 
effect was recently shown to be significant at small film thickness 
by Herrebrugh [5] in a semi-analytical solution for an isoviseous 
and incompressible lubricant. Moreover, Christensen was not able to 
obtain convergent solutions in the final stage of the normal approach 
because of numerical difficulties. 

The present investigation is aimed toward seeking a more effec- 
tive numerical solution for the transient EHD problem which is capable 
of achieving the following: 

1. remove the convergence difficulties at small film thickness, 

2. incorporate the effect of deformation rate, 

3. admit any arbitrary variation of viscosity with pressure, 

4. incorporate the effect of the lubricant compressibility. 

1.2 Previous Investigations 

In spite of the practical significance of the normal approach 
problem, it has received relatively little attention in the literature. 
Before the theories of EHD had been fully developed, Bowden and 
Tabor [6] studied the nature of contact between two colliding solids - 
the collision between a soft metal surface and a steel ball when it 
is dropped from a certain height. Initially, they were concerned with 
the plastic deformation on the dry metal surface by the hard ball 
dropped from a measured height. The initial contact is so small 
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that the impact pressure momentarily reached a value higher than the 
yield stress of the soft metal. The permanent indentation occurred on 
the flat surface when a ball of 1 cm diameter was dropped from a 
height of only 2 cm. To examine the effect of lubricant on the inden- 
tation, they lubricated the flat surface with a viscous fluid, and 
by the electrical conductance method, they detected metallic contact 
and the duration of contact before the ball rebounded. The experi- 
ment with a less viscous lubricant did not give any different results 
compared with those of the dry contact case. The amount of metallic 
contact and the impact time were not altered. However, the experiment 
with a highly viscous fluid showed that during impact metallic 
contact did not occur at all, but the flat surface yielded leaving a 
permanent indentation. This means that the fluid pressure in the 
contact zone at any stage increased beyond the yield stress of the 
soft metal. They explained the phenomenon of surface separation by 
comparing the impact time with the time required to have the fluid 
in the contact region squeezed out completely. If the impact time is 
less than the squeezing time which depends upon fluid viscosity, then 
direct metallic contact is not possible. It is also seen from their 
experiment that the permanent indentation on the lubricated flat 
surface showed a sharp conical shape with the central depth deeper 
than that of the spherical indentation produced by dropping the ball 
on a flat surface from the same height. 

For the first time, Christensen [4] made a theoretical study 
of the normal approach problem of two cylinders in which he considered 
the viscosity of fluid varies exponentially with pressure and the 
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contact surfaces are elastic. He solved simultaneously the two govern- 
ing equations - the transient Reynolds equation and the elasticity 
equation - in time sequence as the gap between the two cylinders de- 
creases. By assuming that the velocity normal to the contacting sur- 
face is uniform within the film, and by employing a direct - iterative 
procedure, he was able to obtain a converging solution for successive 
intervals during the normal approach. However, when the gap becomes 
very thin, the numerical procedure using the direct iteration method 
presents great difficulties and Christensen was not able to obtain the 
convergent solutions in this important region. Moreover, the assumption 
of a uniform velocity is valid only when the film thickness is large 
compared to the deformation. For the small film thichnesses, the local 
normal velocity not only exceeds the center normal velocity but also 
varies drastically along the contact surface. As it will be seen later 
in this work, the local normal velocity at the minimum gap can be order 
of magnitude more than the center velocity. 

Based on his theoretical work, Christensen concluded: 

1. When two elastic cylinders, lubricated with oils whose 
viscosity varies exponentially with pressure, approach 
each other, very high pressures in excess of the maximum 
Hertzian pressure can be developed in the fluid film. 

The elastic deformation forms a pocket shape with the 
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contact. As the film thickness further reduces, the deforma- 
tion tends to flatten out and eventually converges to the 
shape of a Hertzian flat. 

2. For a given load applied to the cylinders, the maximum 

pressure at the contact center depends upon the parameter, 
oE. Harder material and oil with a high a yield a higher 
center pressure during the approach. 

To make qualitative comparisons with his theoretical results, 
Christensen also conducted a series of experiments similar to Bowden 
and Tabor" s work [6] by dropping a ball on a lubricated flat surface 
from a predetermined height. The main objective in his experiment 

was to determine the effects of material properties on the permanent 
indentation on the flat surface. To achieve this, he used several 
pairs of balls and flat surfaces having different material properties. 

He succeeded in proving that under a constant load, the maximum 
transient pressure in the fluid film increases when the parameter aE 
increases. However, he emphasized that this correlation is strictly 
qualitative since the theory is based on the assumptions of an elastic 
cylinder, whereas the actual experiment involves elastic-plastic de- 
formations between a sphere and a flat. 

Recently, Herrebrugh [5], in an attempt of solving the normal 
approach problem of two cylinders, formulated a single governing 
equation by combining the Reynolds equation and the elasticity equa- 
tion. Since he obtained the solution only for the isoviscous case 
which is far removed from the reality of the problem, his solution is 
not complete and his method of solution eventually relies on the 


5 



numerical method, it is hard to see any advantage in his solution 
scheme. The solution of this integral-differential equation - the 
governing equation - is obtained by the method of successive approxi- 
mations with a semi-numerical procedure. He obtained solutions for 
the isoviscous case with the same assumption used by Christensen, that 
is, the normal velocity is uniform within the contact. However, his 
solution only covers regions of high and moderate film thicknesses. 

For extremely thin films, the method of successive approximations 
fails to converge. 

Herrebrugh also noted that as the film becomes small, the ratio 
of the local velocity to the center velocity begins to depart from unity. 
This demonstrates that the assumption of a uniform velocity is no longer 
valid at small film thicknesses. For the isoviscous case at the small 
film thickness where he begins to experience convergence difficulty, 
the ratio of local velocity to center velocity varies from 0.75 to 
1.25. It will be shown in the present work that the pressure-viscosity 
relation has a very strong influence on the ratio of local to center 
velocity at small film thickness. When the effect of variable vis- 
cosity is included in the solution, the local velocity at the edge 
of the contact can be as many as ten times the center velocity. 



CHAPTER 2 - MATHEMATICAL FORMULATION 


2.1 Geometry 

As shown on Fig. 1-1 (a), when the two cylinders approach each 
other along the line connecting their geometrical centers under a 
heavy load, the lubricant between them is pressurized by the squeez- 
ing action of the two cylinders. The contact region where the pres- 
sure is much higher than the ambient pressure is very narrow compared 
with the radius of the cylinder. This fact will be utilized in the 
development of the film thickness formula. The present analysis is 
mainly concerned with the phenomena occurring in this narrow contact 
region during the normal approach of the two cylinders. 

In order to facilitate the mathematical analysis of the problem, 
the contact between the two cylinders as shown on Fig. 1-1 (a) is re- 
placed by the equivalent cylinder with a near-by plane as shown on 
Fig. 1-1 (b). The geometrical requirement for this conversion is that 
at equal value of x the separation between the two cylinders should 
be the same as that between the equivalent cylinder and the flat sur- 
face. 

From Fig. 1-1 (a). 


h = h 1 + R„ 
g o 




2x1/2- 


where h^ is called the geometrical film thickness and h^ is the 
separation on the line of centers. 

Eq. (1) can be expanded to give, 

\ 
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X X 

Since the width of the contact region is very small, (— ) and (— ) 

R 1 R 2 

are both small compared to unity. Thus, by neglecting the terms 
higher than the second power in Eq. (2), we obtain the approximate 
separation between the two cylinders, 


h 


g 




(3) 


Eq. (3) can be rewritten as, 


h 


g 



where 



(4) 


If the radius of the equivalent cylinder is 


R = 


R 1 R 2 

R 1+ r 2 


3 


(5) 


then the geometrical requirement for the conversion from Fig. 1-1 (a) 
to Fig. l-l(b) is satisfied. 


2.2 Governing Equations 

2.2.1 Elasticity Equation 

In the development of the displacement equation a number of 



assumptions can be made based on the relatively small width of 
the contact region where the pressure is higher than the atmos- 
pheric pressure: the contact region is very small compared with 
the radius and the length of the cylinder; the displacement is 
in the state of plane strain; and the tangential displacement 
is neglected because it does not have significant effects on the 
lubricated contact surface. The normal displacement by the pres- 
sure in the fluid film is calculated on the semi-infinite plane 
and then added to the rigid geometrical film thickness. The 
displacement equation is derived in Appendix A and is shown 
below, 

2 o°° 

d(x,t) = - ~ V ^ J P(?> t) in |? - x| d| + C (6) 

— 00 

The constant C is eliminated by including it in the center 
film thickness formula. Due to this constant the displacement 
is not absolute but a relative quantity. 

The film thickness between two cylinders is the sum of the 
rigid geometrical film thickness and the deformations - displace- 
ments - of two cylinders. 

2 2 

h(x,t) = h^(t) + + d 1 (x,t) + d 2 (x,t) + c 1 + c 2 (7) 

From Eq. (7) , 
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h(o, t) = h^(t) 


2(1 - v x 2 ) 


TT E 


1 



p(?,t)MU|d5 


2(1 - v 2 ) 

„ E ~ J p(S,t)jtaU|d5+c +c 2 

2-00 

From Eq. (8) , 

2(1 - ^ 2 ) - 

c l + c 2 = h (o,t) " J p(l,t)i.n| ?|d5 

1 - OD 


( 8 ) 


2(1 - v/) 


rr E 


2 


J P (§,t)Xn|§|d5 


(9) 


Substituting Eq. (9) for in Eq. (7) we obtain 


2 2 2(1 - v n ) 

h(x, t) = h(0,t)+ fgj- + §j£ ^ ^ 


— J p(?,t)jen d? 


2(1 - v 2 2 ) 

rr E 2 



( 10 ) 


Let h(o,t) = h Q (t) 3 which is the center film thickness including 
implicitly the center deformation. 

Define E as, 


l 

E 




(id 


where E^ , and v , 


v 2 are 


Young's modulus and Poisson's 
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Eq. (10) 


ratio of cylinders 1 and 2 , respectively. 

Using Eq. (11) for E and recalling ^ = - 5 — 

J\ 

becomes , 



h(x, t) = h Q (t) + 2 r 


A I 

IT E <). 


p(§,t)Jta 1^1 jp 


d? 


( 12 ) 


which is the film thickness between the equivalent cylinder and 
the flat surface. 


2.2.2 Hydrodynamic Equation 

The inertia force in the flow field between two cylinders 
is negligible compared to the viscous force, which is the funda- 
mental assumption in the derivation of Reynolds equation. In 
the present study the transient, one dimensional Reynolds equa- 
tion is taken as a governing equation for pressure distribution. 
The one dimensional equation is justified by the fact that the 
length of the cylinder can be assumed to be infinite if it is 
compared with the width of the contact region. Further assump- 
tions made in the hydrodynamic equation are: 1 ) the flow is 
isothermal and 2 ) the weight of the cylinder is negligible in 
comparison with the external force. 

The governing equation for pressure distribution is 

fL fph 3 _ 9(Ph) 

Sx Vl2|j, Sx' St 


( 13 ) 


Due to the symmetry of the contact surface at x = 0, the profiles 
of pressure and film thickness are symmetrical at x = 0. 

The boundary conditions for Eq. (13) are 


p = 0 at x = - 00 

^ = 0 at x = 0 


(14) 


Eq. (13) is integrated from x=-xtox=0 using the second 
boundary condition of (14), thus we obtain 



A new variable Q is introduced in order to facilitate the 
use of several viscosity function in the governing equation. Q 
is defined as: 


Q = 1 



(16) 


where p. = — and p, is the ambient viscosity. 
Pi s 

s 

The spatial derivative of Q is 


_d£ _ 1 d(ln Pi) dp 
dx p dp dx 


(17) 


The pressure derivative in Eq. (15) is replaced by Eq. (17), thus 
we obtain 
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( 18 ) 



In the above equation the viscosity term is replaced by ^ 

which is the simple pressure - viscosity coefficient if JJb is an 
exponential function of p. 

Integrating Eq. (18) from x = - 00 to x gives 

Q - - - 12 f [-73 ^ U±> f ^ «] 1- (19) 

- 00 ph ^ -x 

where Q ^ = 0 because at x = - 00 the viscosity is the same as the 
ambient viscosity. 

The value of Q at the film center is 


Q 


o 


12 |l s 




( 20 ) 


The above equation will be used in the calculation of the center 
approach velocity. 

The instantaneous load per unit width of the cylinder is 
the integral of the pressure distribution 


w(t) 


f.CO 

J p(x,t)dx 
- 00 


( 21 ) 


2.2.3 Approaching Velocity 

Since the deformation term in Eq. (12) is the relative defor- 
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mation based on the center deformation which is not known, the 
approach velocity is also the relative velocity, not the absolute 
velocity. However, the relative approaching velocity is incorpo- 
ated in the formation of the present problem because, in general, 
the difference between these two velocities is extremely small in 
the regime of elas tohydrodynamic lubrication. Of course, if one 
would attempt to solve the impact problem of two cylinders like 
the experiment of [6], he should find the absolute velocity which 
plays the important role in the solution of the impact problem. 

Differentiating Eq. (12) with respect to time we obtain 


Bh(x, t) 
9t 


dh o (t) 

at 


4 a r° 
tr e at J 


p(S>t)4n 




d5 


( 22 ) 


It is thus seen that the local approaching velocity consists of 
two terms: the first is the approach velocity of the contact 
center and the second is the velocity due to deformation-deforma- 
tion velocity - which is also dependent on time and varies along 
the contact surface. 

Let v o = -»r and v d ~ A |g J. ^ |? f 1 d? 

then Eq. (22) can be written as: 


ah 

at 


V 

o 


- V . 


(23) 
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2.3 Viscosity and Density Variations 

Both the viscosity and the density of the lubricant are assumed 
to be functions of pressure only. Two types of viscosity functions 
have been used in the present analysis. The first type is the straight 
exponential relation between the viscosity and pressure. This relation 
can be expressed as 

M. = M- e 0 ** (24) 

s 

The second type is the so-called composite-exponential relation between 
the viscosity and pressure. In this relation, the viscosity increases 
exponentially with pressure according to a large exponent in the low 
pressure region and much smaller exponent in the high pressure region. 
Mathematically, it can be expressed as 

p, = for p ^ 

(op + P(p - p )) 

H- = P- g e for p > p 2 (25) 

where p^ = 40,000 psi and p^ = 70,000 psi. 

The viscosity between p^ and p^ is increased asymptotically as 
shown on Fig. 1-2. 

The composite model was first introduced by Allen, Townsend and 
Zaretsky [7]. Their viscosity vs. pressure curve consists of two 
straight lines on the semi-log paper with a discontinuous viscosity 
gradient at p = 55,000 psi. Since this discontinuous gradient is 
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physically inconceivable, before employing their viscosity function in 
the present analysis the discontinuity is removed as mentioned in the 
above paragraph. Their theoretical spinning torque based on this 
empirical equation of viscosity matched excellently with their measured 
torque. The moderation of viscosity increase at high pressure seems 
to be quite reasonable though the exact behavior of the lubricant 
under the dynamic conditions is not known. 

The primary purpose of employing the composite-exponential lubri- 
cant in the present analysis is to understand what effects this lubri- 
cant may exhibit on the pressure, film thickness, load and approach 
velocity. By comparing the two solutions - the one based on the 
straight exponential lubricant and the other on the composite - ex- 
ponential lubricant- one would come up with the plausible conclusion 
on which lubricant model yields the realistic solution in respect to 
pressure and load during the normal approach. 

To find out the effect of O' alone on the pressure profile, the 
two different values of O' in the straight exponential lubricant are 
used in this investigation. 

The density function used in this investigation is 

’ ■ p s( x + r^) < 26 > 

where P g is the ambient density, and and b are the coefficients 
determined from ASME Report C8]. Eq. (26) was originally introduced by 
Dowson and Whitaker [9], 


16 



I 


2.4 Formulation of E las tohydrodynamic Problem 

2.4.1 Coupled Time-Dependent Elas tohydrodynamic Equations 

It has been shown in many previous works on EHD lubrications 
that the solutions for pressure and film thickness must be com- 
patible with each other, i.e., the pressure profile obtained from 
the hydrodynamic equation with a certain film thickness profile 
must be equal to the pressure profile required to deform the 
contact surface to the same film thickness. This demands that the 
hydrodynamic equation and the elasticity equation be solved 
simultaneously at each instantaneous location of the cylinder. 

The two major equations to be solved simultaneously for the 
pressure and film thickness are: 



h = h + 

o 2R tt E 


J P(5,t)jfa ' g - |g| x l d? 


( 12 ) 


2.4.2 Normalization 

Introduce the following non-dimensional variables. 


P = 


. 2_ 


H = £ 


h 

H o = R 2 


, X = ^ , 

3 a 5 


Z = - 


12 ^ s 12 ^ s 

V = v , V = v 

o ER V o V ER 


. D = l 


(27) 
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T 


v 

o 


R 


t. 



9 




9 


w = 



p 





(27) 

cont. 


where a is the Hertzian half-width and the subscript "o’ 1 indicates 
the variables at the film center. 

The normalized governing equations are written as: 


dp 

dX 


/16P V u 
' HZ o p 

' h 3 -p 



dX 


(28) 


H = H + 8P hz 2 X 2 - ( n H ~ ) J P(Z, T) in dZ (29) 

.00 


Eq. (19) and (20) are normalized as follows: 


« - - ( 16 P HZ V o) L [^3 ^ dZ ] “ < 3 °> 

% - ■ ( 16 P HZ V o) r [=£3 ^ «] dX (31) 


The dimensionless load becomes 


W = 



J P(X,T)dX 


(32) 


The dimensionless normal velocity is obtained by differentia- 
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ting Eq. (29) 



P(Z,T) J0n 


Z - X 
Z 


dZ 


(33) 


V = 


- V 

o 



P (Z ,T) Hn 



(34) 


From Eq. (31), we obtain the center normal velocity 



(35) 


2.5 Method of Solution 

2.5.1 Outline of Approach 

Since the pressure and film profiles are symmetrical with 
respect to the center of the contact, it is necessary only to 
obtain solutions for half of a contact. For the present analysis, 
the solutions are obtained in the left half of the contact. This 
half region is further divided into two regions - the inlet and the 
middle region. The division is made in such a way that in the 
middle region the pressure gradient is far steeper than the 
mild pressure increase in the inlet region. 

In the inlet region, the pressure variation is less abrupt, 
and the method of direct iteration can be applied here without 
introducing any convergence difficulties. In the direct itera- 
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tion, the pressure is calculated by the direct integration of the 
hydrodynamic equation for the previously iterated film profile, 
and the succeeding film profile is calculated by integrating the 
elasticity equation according to the newly integrated pressure 
profile. This method is simple and efficient, but is only ef- 
fective for cases of relatively large film thickness. As demon- 
strated by Christensen [ 4 ], for extremely small film thickness, 
the direct iteration fails to yield a convergent solution. 

In the middle region, the system uations are solved by 
Newton- Raphs on method. The solution of the system equations gives 
the pressure correction at every grid point. The Newt on- Raphs on 
method is very effective in solving a system of nonlinear equa- 
tions and usually yields the converged solution in several itera- 
tions. One drawback in the Newton- Raphs on method is the calcula- 
tion of partial derivatives of all the variables in the system 
equations and the inversion of the matrix of which elements con- 
sist of these derivatives. A substantial portion of the calcula- 
ting time for the present problem is expended in the operation 
of the matrix inversion. Details of numerical treatment for the 
inlet as well as for the middle region are given in the next 
sections . 

2.5.2 Integration of Pressure in the Inlet Region 

The integration of pressure in the inlet region is represented 
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by Eq. (30) and is rewritten below: 



( 30 ) 


In the above equation the integral is split into two parts: the 
first in the integral over far left of the inlet region 
(- 00 < X < - X^) and the second is the remaining of the inlet 
region (-X^ < X < - X^) , 

We can approximate the integrand of the first integral, 





d( ft lL. P’) 

BP 


dX (36) 


where we assumed that 


P = 1 , 


H-j. = 1 + 


®hz 2x2 + d ki 


5 (pH) ^ 
5T 


Since the pressure in the inlet region is not high, the normalized 
density is close to unity. is the deformation at X = - X^ 

which is the lower limit of the deformation integral. The defor- 
mation in this region is assumed to be constant. This assumption 
will not produce much error since the approach velocity due to 
the deformation is relatively a small term compared to the other 
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terms in the integrand. 

Regardless of which viscosity model is used in the governing 
equation, the viscosity varies exponentially with pressure in the 
inlet region. Therefore, ^ ~ a • 

Eq. (36) is integrated analytically, 


f^ 1 JL 5(ln H) dy _ 

J 3 3P 


a 
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where 
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The integrated Q written at Kth grid point and time T is 
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(38) 


Once the converged solution for the pressure in the middle region 

is obtained, the integrand in Eq. (38) is assumed to be known 

except density because the pressure distribution in the middle 

region plays the dominant role in determining film thickness 

and approach velocity. In the inlet region the normalized density 

can be approximated to unity for the first iteration. Applying 

the trapezoidal rule for the integration of Eq. (38), we obtain 

Q__ • Then P Tr in the inlet is determined from Q TT as : 

x K,m K,m x K,m 
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Q K,m 1 " n, 


K,m 

-op. 


= 1 - e 


K,m 


(39) 


Thus the pressure equation in the inlet region is 

P* = - 7 : An(l - Q„ m ) (40) 

Kjin Ctf K,m 

2.5,3 Calculation of Deformation 

The deformation for an arbitrary pressure distribution can 
not be determined by the straightforward numerical integration 
because the integrand in the deformation equation becomes singular 
at X = Z, Care must be exercised in the formulation of the nu- 
merical integral formula by which the singularity at X = Z can 
be removed. 

The detailed derivation of the quadrature formula for the 
singular integral kernel is presented in Appendix A and the 
quadrature formula is written below, 

K0-2 

r - 1 tvi (-v zj ) 

"TCI j=l,3 ,5 


P. K„ 
j+l,m 2 


(- X K,- Z j) +P J+2,m K 3(- X K ,- Z 3 ) - 


(41) 


where 
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(42) 
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and 
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(45) 



(46) 

(47) 

(48) 


Since the pressure profile is symmetrical about X = 0, the second 

half of the deformation integral can be approximated in the same 

form of Eq. (41) by changing -Z. to Z. in X and K q , thus 

J J -t » * 3 


KO-2 


f\<Z>fa|z - Xjdz = l + h(-\ Z S ) ] 

“TCI j=l,3,5 
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and following the above procedure we obtain 


X KI Ko_2 

J P m (Z)to|z|dZ -l + ] 

"nil j=i,3,5 


+ *i+l,l*2(ho,- Z i) + ^ho, Z S ) J 
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where = 0. 

For the convenience of differentiating m with respect to 

P. , K- , K n and K 0 are rearranged in such way that P. has a 
j ,tn 1 2 3 ° y J,m 

single coefficient RC-X^ * Z^): 



j = 1 
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where 
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The final form of the deformation equation is 
Ko 

D K,m-' C 3 l E (-V Z j) P 3.m 
3 = 1 , 2 ,- 


16P, 
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( 52 ) 


(53) 


2.3.4 Elastohydrodynamic Equation in the Middle Region 
Eq. (28) written at Kth grid point and time T m is 



(54) 


The derivative 



terms and can be approximated by 
quadrature as 


(54) may be split into three 
the Lagrangian three point 
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and 


H = H + 8P„_ 2 X_. 2 
g K>m o,m HZ K 


(59) 


The first two terms on the right hand side of Eq. (55) can be 

grouped together and expressed by Y m (-X K ) in which all the 

variables were determined in the previous time steps. Therefore, 

V (~X^_) is not a function of P . 
m K j ,m 

After rearranging the integrand in Eq. (55) to a pressure 
dependent term and a pressure independent term, Eq. (55) may be 
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written as 
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where A X. = X.^- - X. 1 K + 1 ^ i ^ Ko-1 

i 1+1 1-1 

= X_- - X. i = Kj Ko 

l+l l 


Substituting Eq. (53) for D_^ ^ in Eq. (61) and rearranging 
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where 
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The integral term and the deformation terms in Eq. (54) are 
replaced by Eqs. (62) and (52), respectively. The discretized 
form for Eq. (54) at -X^^^ can t* 1118 written as 
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(63) 


Eq. (63) is one of the typical equations in the system equations. 

If ¥ (P) is written at every mid point between grid spacings in 

the middle region, there are N equations with N unknown, p , 

K,m 

where N is the number of grid points in the middle region. 

Applying the Newton-Raphson technique to the system equations, 
we obtain 
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(64) 


{yp)} 0 ?^ 


(n+l) 

i i 


m* 


L 1 ^] 



where ^ j* and J represent a column matrix and an N x N matrix, 

respectively, and A* indicates partial derivative is to be taken 

with respect to P . n is the level of iteration, 
m 

From Eq. (64) we obtain 
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The right hand side of Eq. (65) is assumed to be known from 
the lower level iteration, and "|a P is defined as 


{a P V n+1 ) = {p |( n+1 ) _ -jp |( n > 

m> ^ mJ 
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The elements of the matrixes in Eq. (65) are detailed in Appendix 

B. 

The center approach velocity and the load at time T are 
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where 
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for the straight exponential lubricant 


and 


-fSP s + B(l-P s )} 

Q q = 1 - e for the composite-exponential 

lubricant. 


The film thickness written at Kth grid point and time T m is 


1 + 8p hz \ + d k. 
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(69) 
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2.5.5 Outline of Numerical Procedure 

For the computational convenience, it is assumed that the 
center pressure is constant while the value of load varies as the 
cylinder approaches the flat surface from a high point. The 
calculations are performed to obtain the several series of the 
solutions in which each series represent the solutions at various 
center film thickness with a fixed center pressure. 

The best approach to the problem is to obtain analytically 
the pressure distribution for a high center film thickness by 
neglecting the deformation term in the hydrodynamic equation, and 
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at each time step the center film thickness is reduced a certain 
amount and is kept constant. 

Written below are the precedures of numerical calculation at 
each time step: 

1) At the first time step analytically obtained pressure 
distribution is used as an initial guessed pressure. 

From the second time on, the initial guessed pressure 
is determined by linearly extrapolating the previous 
pressure distributions. 

2) Using the initially guessed pressure distribution, the 
film thickness, density and viscosity are calculated. 

Then the approach velocity is determined from these 
values. We set up system equations (63) to obtain the 
pressure correction terms in the middle region. Once 
the pressure distribution in the middle region is 
corrected by ^A P J - , the inlet pressure profile is de- 
termined by linear interpolation with the factor 

P A ( n+ l)/p ( n ) w h ere p ( n+ l) £ S obtained from 

KA,m KA,m KA,m 

the system equation. The film thickness is calculated 
using the newly obtained pressure. 

3) If the converged solution for the pressure in the middle 
region is obtained, Eq. (38) is solved for the inlet 
pressure and the center approach velocity ^ is de- 
termined by Eq. (67). Now the overall pressure distribu- 
tion is checked for convergence. If it has converged, the 
load is calculated by Eq. (68) and one moves to the 
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next time step. Otherwise, the above procedures (2) 
and (3) are repeated until the converged solution is 
obtained. 
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CHAPTER 3 - DISCUSSION OF RESULTS 


3.1 Introduction 

The results of the present study are presented as a series of 
curves for pressure, film thickness, load and approach velocity cal- 
culated at a prescribed center pressure and at successive reductions 
of the center film thickness. 

The pressure and film profiles for various parameters at successive 
stages during a normal approach process are plotted for the left half 
of the contact region. The integrated load and the approach velocity 
during each normal approach are plotted against the center film thick- 
ness or the minimum film thickness. 

3.2 Pressure Profiles 

Shown on Fig. 1-3 to 1-13 are the series of the pressure profiles. 
Each figure displays the change in pressure with film thickness as 
the cylinder approaches the flat surface for a given center pressure. 

The range of the center pressures employed in the present study is 
from 2.5 X 10 4 psi (1.723 X 10 8 N/m 2 ) to 1.5 X 10 5 psi (1.034 X 10 9 N/m 2 ) 
which are typical maximum stresses encountered in concentrated 

contacts . 

In general, the trend of change in pressure with respect to the 
center film thickness is qualitatively similar for all cases, namely, 
at high film thickness the pressure level decreases steadily through- 
out the contact region with decreasing film thickness until it reaches 
a stage when the integrated load becomes a minimum. After this 
stage the pressure in the middle region reverses its trend and begins 
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to rise, but the pressure in the inlet region still continuously de- 
creases as the center film further decreases. In all cases, the pres- 
sure rise is confined within a small fraction of the Hertzian half- 
width, and it does not appear to reach the Hertzian semi-elliptical 
shape. 

For the straight-exponential lubricant, the pressure-viscosity 
coefficient, a 9 has a marked influence upon the pressure gradient near 
the center of the contact. For example, Fig. 1-9 shows that the pres- 
sure gradient for a = 12.8 at the center is far steeper than that 
appearing in Fig. 1-5 for a = 9.5. 

The change in the center pressure also produces a very strong 

effect upon the pressure gradient at the center. A higher center 

pressure produces a sharper pressure spike at the center. The effect 

becomes increasingly stronger at higher center pressures. For example, 

9 2 

at center pressure equal to 150,000 psi (1.034 x 10 N/m ), the pres- 
sure gradient gradually tends to become infinite. The existence of 
such sharp pressure spikes in practice appears to be highly question- 
able, since the shear stress would also become incredibly large under 
these circumstances. It appears very unlikely that the fluid can 
withstand such high shear stresses, particularly in the light of 
recent work on traction studies [lO], [ll], and [l2] which demonstrate 
the existence of a limiting shear stress for any lubricant. In the 
vicinity of this limiting shear stress, the fluid behaves in a non- 
Newtonian fashion, and an increase in shear rate has little effect on 
the shear stress. 

The effect of the non-Newtonian behavior can be accounted for indirect- 
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ly by introducing the so-called composite-exponential model for the 
lubricant viscosity. This was demonstrated by Allen et al [7] in a 
spinning torque study. The resulting pressure profiles using a com- 
posite-exponential model similar to that in [ 7 ] are shown in Fig, 

1-10 to 1-13. These curves show considerably different features com- 
pared to the pressure curves for a straight exponential lubricant. 

For example, the pressure gradient is much more moderate near the 
contact center, showing the absense of a pressure spike which is so 
characteristic for the straight exponential lubricant. Moreover, the 
steepness of the pressure gradient near the contact center is not 
influenced greatly by the increase in the center pressure. For example, 
there is very little difference in the pressure gradient between 
Fig. 1-10 and Fig. 1-13 at the same film thickness. 

It should be emphasized that the results for the composite-expo- 
nential lubricant are intended to show the qualitative effect of the 
reduction of pressure-viscosity coefficient on the characteristics of 
pressure and film profiles. These results should not be used quanti- 
tatively for design purposes. 

3.3 Film Thickness 

The film thickness profiles are plotted in conjunction with the 
corresponding pressure profiles in Fig. 1-3 to 1-13. At the early 
stage of normal approach, a pocket is formed elastically at the contact 
center, and its shape does not change much for subsequent reductions of 
the center film thickness. The pocket depth defined as the difference be- 
tween the center film thickness H and the minimum film thickness, is depen- 
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dent upon the center pressure for a given lubricant. A higher center 
pressure produces a deeper pocket. 

When the center film thickness decreases to a certain level, a 
quite different phenomenon occurs. At this point, the normal approach 
velocity at the center suddenly drops almost to zero, while 
the local approach velocity elsewhere in the contact continues. 

This condition produces a deeper pocket during the final stages of the 
normal approach. In all cases investigated, the growth of the pocket 
persists all the way down to the very end when the edge of the contact 
at the minimum film thickness point practically touches the opposing 
surface. For perfectly smooth surfaces, the point of the minimum film 
would eventually form a seal and the lubricant inside this point 
would be trapped. Thus, by including the local approach velocity 
in the analysis, one can show that both the pressure and film thick- 
ness profiles never reach the semi-elliptical Hertzian shape as sug- 
gested by Christensen in [4l. Instead, the pressure remains to be 
confined in the center region, and the surface deformed into a pocket 
inside which a portion of the lubricant is entrapped. As shown in 
these deformation shapes, the center pressure has a definite influence 
upon the depth as well as the width of the pocket. In general, the 
pocket becomes deeper and wider as the center pressure increases. 

The pocket formation is more pronounced for the case of the com- 
posite exponential lubricant. The pocket depth is somewhat greater 
than the corresponding case for the straight exponential lubricant. 

The change of the pocket shape during normal approach is qualitatively 
similar to that for the straight exponential lubricant. At the last 
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time step when the minimum film thickness H, is less than 5 x 10 

the pocket depth increases rapidly while the location of the minimum 

film thickness moves slightly toward the outer edge of the contact 

region. The highest value of pocket depth for all cases investigated 

5 9 2 

occurs at a center pressure, P q = 1,5 x 10 psi (1.034 x 10 N/m ), 
with the composite exponential lubricant. The value of the maximum 
depth exceeds 30 x 10 and there is practically no significant 

pressurization outside of the pocket. It is thus expected that during 
the normal approach of two cylinders the pressurization is effectively 
contained inside the pocket and that the width of the pocket is approxi- 
mately one-half of the Hertzian contact width based on the same center 
pressure. 

3.4 Load 

Shown on Fig. 1-14 are the load vs. center film thickness curves 

at a constant center pressure for the straight-exponential lubricant. 

In general, the dependence of load on the pressure-viscosity coefficient 

a and the center pressure in the present analysis confirms Christensen's 

conclusions: first, for a given center pressure, the load is strongly 

dependent upon the pressure viscosity coefficient, i.e., the higher 

a produces much smaller load* For example, the load for a - 12.8 

8 2 

and = 100,000 psi (6.894 x 10 N/m ) is approximately equal to the 

8 2 

load for a * 9.8 and P q = 25,000 psi (1.723 x 10 N/m ); and second, 

once the center pressure is sufficiently high, the increase in load 

is negligibly small for further increase in center pressure, i.e., 

the load becomes insensitive to the center pressure. As described 
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before in Section 3.2, this insensitivity of load to the increase in 
center pressure is caused by a strong pressure-viscosity coefficient 
O'. Thus, one would expect that if the increase in viscosity with pres- 
sure is milder, the load becomes more dependent upon the center pres- 
sure, as will be seen in the results of the composite-exponential 
lubricant. 

Also in Fig. 1-15, a quantitative comparison is made between 
the load curves obtained by Christensen [4] and those calculated from 
the present analysis. On the right side of the minimum load, the two 
theories shows fairly close agreement, the present analysis yielding 
a slightly higher load. This slight discrepancy in load is attributable 
to two effects: first, the approach velocity in the present analysis 
is higher than that in [4] where the local deformation velocity is 
neglected, resulting in stronger squeezing action on the fluid by the 
cylinder, and second, the effect of the compressibility of the lubricant, 
which was also neglected in [4]. On the left side of the minimum load, 
the effect of the local deformation velocity becomes very important, 
and the present theory gives considerably higher load than Christensen's 
results. Furthermore, there is also considerable difference in slope 
between the two results. The present theory predicts a much steeper 
slope on the left side of the minimum load, indicating that there is 
virtually no reduction in the center film thickness while the minimum 
film thickness steadily drops to zero as shown on Fig. 1-15. 

It should be noted that the maximum load obtained in the present 
analysis is substantially less than the corresponding Hertzian load 
based on the same center pressure. This result directly contradicts 
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Christensen’s conclusion that the load increases to the Hertzian load 
as the minimum film thickness decreases to zero. 

As shown on Fig. 1-18, one may find the variation of center pressure 
at a constant load during the normal approach of the two cylinders from 
Fig. 1-15 and 1-17. If a horizontal straight line is drawn at any 
specific load on Fig. 1-15 or Fig. 1-17, depending upon the lubricant 
used, the change in P q with decreasing center film thickness can be 
determined from the intersection of the straight line and load curve. 

The center pressure gradually increases with decreasing center film 
thickness, and then increases abruptly to the maximum value; the maxi- 
mum is much larger than the initial p^. The center pressure finally 
decreases rapidly for further decrease in center film thickness. 

In Figs. 16 and 17, results of the composite-exponential lubricant 
show that in general, the loads are much larger than the corresponding 
loads for the straight-exponential lubricant. The change in load with 
the center film thickness, or with the minimum film thickness, is some- 
what moderate. No abrupt increase in load is seen. The most noticeable 
effect produced by the composite-exponential lubricant is the relation- 
ship between load and center film thickness. The load is strongly de- 
pendent upon the center pressure. 

3.5 Approaching Velocity 

As mentioned in Section 2.2.3, the center approach velocities 
shown on Fig. 1-19 are not the absolute velocities - the 
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velocities of the approaching cylinder center they are the relative 
center approach velocities, i.e., the time derivative of the 
center film thickness. However, it is known that in the normal approach 
problem of EHD lubrication the difference between them is negligibly 
small. 

It is apparent from Fig. 1-19 that the center approach velocity 

V decreases with decreasing center film thickness at a constant center 
o 

pressure, and the rate of reduction in V is a function of H and P . 
r ’ o o o 

In the region of high H q , the center approach velocity approximately 

varies with the square of the center film thickness for a given center 

pressure. This trend agrees with that predicted by the normal approach 

solution between two rigid cylinders. This parabolic relation between 

H and V ceases to exist as H is reduced to a certain value depending 
o o o * ° 

5 8 2 

upon P q . For example, for P^ = 1.25 x 10 psi (8.617 x 10 N/m ) and 
H q approaching 3 x 10 V q decreases rapidly for further 
decrease in H q . For low center pressure, this transition occurs at a 
much smaller value of H q . The rapid reduction of the center 
approach velocity for high center pressure can be explained by con- 
sidering the flow quantity through the gap between the bump and the 
flat surface* The gap is not more than 10 microinches so that the 
lubricant flow through this gap is very small; consequently very little 
squeezing on the lubricant is necessary to maintain a constant P q . 

It is interesting to note that the center velocity V q required to 
produce a high center pressure P q at a constant center film thickness 
H q is considerably lower than that for a lower P^. This trend directly 
opposes that based on the rigid cylinder theory for which a greater 
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P q requires a high center velocity V q at a same center film thickness 
H q . This discrepancy can be accounted for by the deformation effect. 
At a higher pressure, the contact region is larger, the squeezing 
action is thus much more effective; and it requires a smaller center 
velocity to produce the required center pressure. 

Fig. 1-20 shows the ratio of local approach velocity to center 
approach velocity vs. H/W for three points of the contact region 
X = -0.25, -0.5 and -0.75. For the sake of comparison, typical data 
from [5] are also shown on Fig, 1-20, As expounded in Section 2.2.3, 
it is known that local approach velocity varies along the contact 
surface and the most severe variation occurs when the film thickness 
is very small. The data from C5 D is based on the assumption of iso- 
viscous lubricant, which shows the variation of local velocity is 
relatively small compared with that for the lubricant of variable 
viscosity. This comparison clearly indicates that it is much more 
difficult, sometimes almost impossible, to obtain the converged solu- 
tion when the center film thickness is small because controlling the 
local velocity numerically between two successive iterations is very 
difficult. 
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CHAPTER 4 - SUMMARY OF RESULTS 

It has been found that the full solution of the normal approach 
problem of two elastic cylinders, with a compressible lubricant between 
them whose viscosity varies exponentially with pressure, can be obtained 
by solving numerically the coupled transient Reynolds equation and the 
elasticity equation using a combination of direct iteration and Newton- 
Raphson method. 

The results show that: 

1) In general, the pressure profile for the straight exponential 
lubricant shows a sharp spike near the contact center; a 
higher center pressure or a higher pressure-viscosity coef- 
ficient results in a steeper pressure profile at the contact 
center. However, for the case of the composite-exponential 
lubricant the steepness of the pressure profile at the con- 
tact center does not depend so strongly upon the center pressure. 

2) For all cases studied, a pocket is formed elastically on 
the cylinder surface near the contact center during the 
early stage of the normal approach, and it remains without 
much change in its shape until the final stages of the normal 
approach, resulting in a quantity of lubricant inside the 
pocket being entrapped. Thus, the film profile never reaches 
the semi-elliptical Hertzian shape as suggested by 
Christensen [4], The depth of the pocket is dependent upon 
the center pressure for all cases investigated. In compari- 
son, the pocket depth for the compos ite- exponential lubri- 
cant is much deeper than the corresponding one for the 
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s traight-exponential lubricant 


3) In general, the load increases very rapidly from its minimum 
value with virtually no reduction in the center film thick- 
ness. This result can be attributed to the fact that the 
entrapped lubricant inside the pocket is effectively pres- 
surized further by closing the gap between the minimum film 
thickness and the flat surface. This pressurization, in 
turn, deepens the pocket depth further. Thus, for all 
cases investigated, the load never increases to the Hertzian 
load based on the same center pressure as the minimum film 
thickness decreases to zero. In contrast to the cases for 
the straight exponential lubricant where for a sufficiently 
high center pressure and at any given center film thickness 
the load is insensitive to the center pressure, the load fo 
the composite-exponential lubricant is strongly dependent 
upon the center pressure. 

4) At early stages of the normal approach, the local approach 
velocity does not deviate from the center approach velocity. 
However, during the final stages, the ratio of local velocity 
to center velocity greatly exceeds unity, indicating that the 
center film thickness is almost constant while the film 
elsewhere continuously decreases. For a given center film 
thickness, the center approach velocity required to produce 

a higher center pressure is considerably lower than that 
for a lower pressure. This trend is more pronounced at the 
final stages of the normal approach when the deformation 
overtakes the geometrical film thickness. 
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APPENDIX A 


QUADRATURE FOR INTEGRATION OF ELASTICITY EQUATION* 

Referring to [13 ] for detailed derivation, the normal displace- 
ment for any x on the surface of semi-infinite solid due to vertical 
forces is given by 

p ^ Iz - X I 

V x > - J V z) * n - L firr L dz <Aa) 

" Kl 

where the symbol |z - x| represents the positive distance between the 
force element at Z and the point of interest at X as shown on Fig. A-l. 



Since the integrand is singular at X = Z, the numerical quadrature 
formula should be developed in such a way that the singularity at X = Z 
can be removed. It consists of approximating the function P by a para- 
bolic polynomial in each subinterval, performing the integration in 
closed form in the subinterval, and summing over the whole region of 
integration. 
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We subdivide the right half of the contact region into N sub- 
intervals, requiring that the widths of two consecutive subintervals 


equal and assuming the pressure distribution is known. Then 


I(X) = 


J P m (z)fa|z - xj 


N 

dZ=V j ( X ) 

J-l 


(A. 2) 


where 


p J j+l 

I (X) = J P m (Z)^n|z - xJdZ 

J *7 


(A. 3) 


Z. 

3 


The parabolic representative of the pressure distribution in the 

subinterval [z., Z is 
J J +1 


P m< Z > - [ ( Z - Z j+ l/ 2 ) ( Z - Z J+l) P jjm - 2 ( Z - Z j) ( Z - Z j+l) Vl/2, 


m 


+ ( z - z i) ( z - z J+ i/ 2 ) Vi, j /2 V 


(A. 4) 


where 




Z. 

J 


) 


From (A. 4), 
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• r ;( z j) - (- 3r J>n + 4 f j + 1 / 2 ,„- 

P j + l,m - r ;( Z j«) - ( p jjn - 4P j+ l/2,m + VlJ' 21 ), 

P j,m ” P j+l,m " ( p j,nf 2P j+l/2,m + p j+l,m/ /A j 2 


We integrate Eq. (A. 3) by parts several times to obtain 


V X > = [ p B »)0 £ _1 {*>U - *l} - P;(Z)D f " 2 {*.|z - X|} 


z 


P m (Z)D f " 3 { Xn l Z - X l} ^ 


3+1 

3 


where 


D f _1 {in|z - X|} = J J8n(z - X jdZ 
D f " 2 {£n|z - x|} - JJ* jGa|z - x| dZ 
D f " 3 {^n]z - X |} = JJJ Anlz - x|dZ 


Thus 

B f ~ l {jLn\z - X lj = (Z -X) hx ]z - X I - (Z - X) 
D f " 2 {in|z - X(j = | (Z - X) 2 Xn | Z - x| - | (Z - X) 2 


(A. 5) 


(A. 6) 


(A. 7) 


(A. 8) 
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D f ' 3 {jJn|z - x|} = (Z - X) 3 Xn|z - x( - (Z - X) 3 
Substituting (A. 8) in (A. 6) and some manipulation yields 


(A. 8) 
cont. 


IjOO = [(z - X){j2n|z - x| - l} {p m (Z) - \ (Z - X)P^(Z) + (Z - X)V(Z)/6} 

Z x 

+ (z - X) 2 P*(Z) - 5(Z - X)P"(Z)/9}/4~! (A. 9) 


z. 

J 


Le t Uj - Zj - X K> u j+1 - Z J+1 - X and Sj - f - f j 

with these variables and noting that at the end points of each sub- 
interval in the interior of [- X 0], there is exact cancellation of 

Kl, 

the P^(Z) contribution, Eq. (A. 9) is rewritten as: 

I.(X^) = (p. S.- P' S...J - \ P'fu. (s. - \ U. 2 ) 

J K ' J,m J J+l,m j+1' 3muj\j 6 j / 
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" u j+l ( s j+l ■ 6 u i+l ) ] 


j+1 


(A. 10) 


Substituting Eq. (A. 5) for P 1 and P n in (A. 10) and summing over 
the entire interval. We obtain, 

r X KI «- 2 r 

J P m C2>-^|z - ^IdZ - ^ 


1,3,5 


+ V,A<w! 


(A. 11) 


where 
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K i(Y Z j) ’ 26 - ^ 3S j + S j+ 2 ^ ~ OR 2 ~ " 1 ] ’ 

3 j 

K 2(\ Z j^h( S j +S j + 2) + ^2 > 

j __ 

^^Sc, 2 ^ = 2®T ^'’ S j " 3S j+ 2 ^ " TT* 2 + U j+2[ An ^ U j+ 2 ^ " 1 ]’ 


36. 
J 


and 


S. - u.(s. - ~ u . 2 ) - u (s._,„ - T U 2 ) . 

J J'J 6 j / j+ 2 \ J +2 6 j +2 / 


The quadrature formulation for the singular kernel in the integrand 
written here is exactly the same as that of Ref. [14] . 
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APPENDIX B 


CALCULATION OF MATRIX ELEMENTS IN EQ. (64) 

For (convenience, Eq. (63) and (64) are rewritten below 


BP 


R+l/2 ,m 


^m^K+1/2 ) 


• j- / , xii 

f P K+1 , m" P K , m ^ ( 1 + A l P K+l/2 ,tiA 

= ' ^ } ( -<* w;y } 


KO 

(h - c~ V r(- x_. 1/9 - Z.) P. [ 

v g K+l/2,m V ^ +1/2 > J ' 


j=l 


KO 


- (8P„_V ) { y I v (-x.) - 0 ) p. (h - l) "I AX. 

' HZ o,nK <■ L> L m' i m i,m' g. m * J l 

i=K+l/2 


'i,m 


u) C- 
m 3 


KO 

I 

j=l 


K.“ 



(63) 


U (p)} + {ap j- Fa*? 

* m J v. mJ L m 



(64) 


The calculation of the matrix elements in [A-Y (P) ] involves the 

m 

differentiation of {Y (P) } with respect to {p }. Before differentia- 

m m 

tion, Eq. (63) is rewritten in the following form: 


vw4 (w \J 

K ^K+l/2 ,m 


( 4P HZ V o,m/' (' 


I K+l,m + I K,m 


) 


(B.l) 


k 
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where 


KO 


KO 


L, = Y [y (-x.) - cu p. (h - l)l AX.+ .C.T l(-X -Z.) P. 

K Zj L m 1 m x,m \ g. SJ 1 m 3 Lj \ k 3 j / j ,m 

i=K 1 » m '-•« 


j=l 


(B.2) 


and 


The variables, m ant * *k+ 1/2 m’ are ex P resset * as th e average 

of the two values at -X^ and as : 


I K+l/2,m 2 ( I K+l,m + 

*St+l/2 ,m = 2 ( H K+l,m + ^m) 


(B.3) 


(B.4) 


The ^K+l/2 m anc * ^K+l/2 m are ta ^ en as a f unct i° n of the average 
pressure, (P K+ ^ m + m ) and expressed below: 


!*(> 


2 " V P K+l,m + P K,m 


) 


pK+1/2 - m ’ uhL + p v ) 

2 1\ K+l,m K,m' 


(B.5) 


“*.+1/2 ,m * exp B " ( r K+l,m + \,m) ] 


(B.6) 


The derivative of the variables in Y (P Tri1/0 ) are derived below: 

m K+l// 


sr*® - Vs 1 '-*!,- V - (I) a*K, n - » 

J K,m 


(B. 7) 


where 
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for j 2: K 


8=1 

s 


6 =0 for j < K 

s J 


and 


K0 


if j / KA, L(-X^ - Z.) . V p. R(-X. - Z .) AX. 

K, j Aj x,m l, J x 


i=K 


K0 KA 


if j = KA, L(-X^ a - Z.) = ) Y p. R(-X. - Z.) C 1?m -) AX. 

TKA, j' Li i-> x,m x, y ^V A J i 

i=KA j=l ’ 


In this way, we can take into account the effect of the pressure 

distribution in the inlet region on D__. - the deformation at the 

dividing point between the inlet and middle region, since D„. is 

KA,m 

strongly dependent upon the inlet pressure distribution. 

If j = K or K + 1, then 


( e „n 3 f P + P y 

3P . 3P . V exp L2 \ K+l ,ni P K,raJ / 


9p. 


2 “ M 'K+l/2,m 

K-KL/2 ,m = 5 f 2 B ( P K-fl,m + P K,m) 1 

3p j,m ap 3,m 1 + \ A l( P K+l,m + P K,m) 


B 


2 + A 1 1 P„. . + P_ ) 

1 \ K+l ,m K,m/ 


(B.8) 


CB.9) 


Since the deformation, D depend upon the overall pressure dis- 

K,m, 

tribution, the derivative of D„ with respect to any P. exists. 

K,m v J j,m 
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Thus 


'" IT- 2 -' 3 * - ( 3 VYj) (^fl,n + O 

J > m 


(B. 10) 


where 


“k 


,J - R (- Xk+ 1,- Z j) + E (- X k,- Z j) for j f* KA, KA + 1 


KA 


“k,3 - I { Z i) + E ('V Z i) ] (p^f ) } for 3 - “• “ + 1 

KA,m 


i=l 


The reason for summing the products of the deformation kernel and 

the inlet pressure ratio over the entire grid in the inlet 

region is to take into account the effect of the inlet pressure 

distribution on D TrA + D TrAr1 

KA,m KA-fl,in 

Using Eqs. (B.7), (B.8), (B.9) and (B.10), the derivative of 

^m^K+1/2^ is written as: 


9Y (P 
m 


dP 


r^M! s) kj (w ( - K X~ k 




- ( 4 E HZ V o,m) { [Vst 1 !- V Z j) + L (‘ W Z j) ] - ( 6 


-) 
s mS 


^2+A (p “ +P K / ^ K -”" ^ ^ + ( ' Hk+1 ’»' ^ ^ } 

l x K+l,m K,m/ 

+ (6 ) { L . P ) I" a** 

l V 8A5L. I \ K+l ,m *K,mJ L 


'2 + A 


l^ P K+l,m + P K,m) 


(B. 11) 
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+ & a )©]-fe a )( 0 } 


K+l/2,m 


^+1/2, m' ' S< 


(B. 11) 
cont. 


where 


8=0 

u 


6 = 1 
u 


j / K, K+l, 
j 7* K, K+l, 


and 


6=1 

g 


j - K+l, 


8 = -1 
g 


j = K . 


Eq. (B.ll) is one of the typical matrix elements. The expanded 


form of Eq, 

(64) is 


»V P KA> 

M (P„ A ) 
m KA 

(*V A ) 

m KA. 

BP KA,m 

^ P KA+l ,m 

9P K0-l,m 

3 VW 

SI .V> 

*VW 

KA 3 m 

Bp 5 

r KA+l,m 

dP K0-l 3 m 

i 

t 

li 

i 

l 

s 

1 

l 

i 

i 

9 V P K0-1> 

3 V p ko-i> 

8 VW 

SP KA,m 

^KA-fl,m 

dP ’ 

K0-l,m 


AP. 


KA,m 


AP 


KA+1 ,m 

i 


AP. 


K0-1 ,m 




VW 




Y n/ P K0-1^ 


(B. 12) 
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The pressure correction term at the contact center, AP Trrk , is 
r KO ,m 

not necessary since the center pressure is assumed to be constant. 
The center velocity V q ^ is kept constant during the calculation of 
the pressure correction terms. The center velocity is recalculated 
after the converged solution for the pressure distribution in the 
middle region is obtained. 
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appendix c 


COMPUTER PROGRAM FLOW DIAGRAM AND FORTRAN LISTINGS 
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PROGRAM ElASTO (INPuT, OuTPjT, PoNCH, I APE5= I NPJT , TAPEb=OuTPbT ) 
COMMON / A/ KF* AlPH, EN, ED» HE, HE 1 , ALPH1, MM 

COMMON iv I » KA* KKl* KKF , KF2, KF3, V , KKE* lb* N» NN» PO* K* .-i * 

1 P t*i * PC » Kb , EP , EU , w , W ft » .V ri , RA * t\U * PNK * PD I » Xh , OX » u I , 0 , 
2 HG»F(*P*VIo»VISD»DcN»DlND» T ,OT *Vb,DP»HO»bE.»A*L5»C*SL , »Ai - i*Pb* 

30MEG3 » OMEG2 »OMEG 1 » DENT ,DD*bb*C4*C3 , I 1 *w » P A I » BE T » P 1 »P2»VIbl * 

4V I S2 » V I o3 » DP V » C6 » EBP » bUB V » CC » c c. » P KA * DPR , K 
D I MENS ION XH (31) ,DX(31),0(Pl»3),HG(3i,3) ,H(pi,3),0I(bl), 

IP ( 51 » 3), VIS ( 31 ) * VI bO (31) ,DuN (51,3) , DEND (51»3)»T(35)*V(35), 

2DT ( 35 ) , VD( 51 ) *DP ( 51 > »HO( 50 ) *.-. < 35 ) *SE ( 51 ) ,A( 31 ) ,d( 51 ) *CC( 51 ) , 

3 F ( 5 1 ) ,R(51,51 ) » BE T (51,51 ) *60(51') » PS ( 5 1 ) » OMc.b 3 ( 3 ) » UMEb2 ( 3 ) , 

40 MEG 1 ( 3 ) , DENT (51), 0(51, 51 ) » PA 1(51) »SJBV (51), C(2, 3, 51), EE (51) 

D I MENS ION TEM ( 5 1 ) , N I ( 3 1 , 4 ) 

R E A D ( 5 » 2 u 0 } KI,. KF, MI, MF , ED, EN ,' ALPH, PH 
READ(5,203) EP, ED, HOI, RA 
READ (5,206) ( X H ( K ) , < = K I * K F ) 

READ(5,2d5) ( ( Q ( K , J ) , J = KI,kF), n = K 1 , KF ) 

KKl=KI+l 
KF2=KF— 2 
No = 3 
KF3=1 
KkF=KF- 1 
E = 3 3 • 6E+6 
H01=20o • Ob— 6 

EH = FN*1.25iED = EU*1.23bALPH = ALPH*1.25j.PH=123000.U 
WRITE (6,2 17) ED, EN , ALPH, PH 
WR I TE ( 6 ,841 ) ( K,XH( K ) ,K = K I ,KF ) 

N r = 2 S N I = 3 
KA = 33 
PnZ=PH/E 

C4=16.u*PHZ b Cp=C 4*PHZ b Cb = C 3 / 3 • 14 
WRITE ( 6 , 2 ■ j 5 ) PHZ , C4, C5 , C6 , E 
IP(NT.EQ.l) l7, 39 
27 READ (3,^03) ( P ( < , 2 ) » K = K I , KF ) 

READl5,2u5) (P(K,1), K=KI, KF ) 

READ ( 3,231) H ( KA , 1 ) , T ( 1 ) 

READ ( 5,231) H ( KA * 2 ) , T ( 2 > 

I I =2 b Ib = 2 b MM=1 4> M=1 a> CaLl DVD b Call hDT 
MM=1 $ I S = 2 b M = 2 b CALL OVD b CALL HDT 
W R I T E ( 6 , 6 6 8 ) (K, P ( K , 2 ) , K = K I , KF) 

WRITE (6, 2-6) H ( KA , 1 ) , T(l) 

WRITE(6,2o6) H ( KA , 2 ) , T ( 2 ) 

WR ITt(6»668) (K, H(K,1), K = KI, KF ) 

w'R I TE ( 6 ,668 ) (K, H(K,2), K = Rl, KF ) 

wKITE(6»668) (K, D ( K , 1 ) , K = K 1 , KF) 

WR I TE ( 6 ,668 ) (K, D ( K , 2 ) , K = KI» KF ) 

CALL INTEG1 ( XH ( Kl ) >XH ( KF ) ,2 ,P ( 1,2) , KF , VA L J t , I ER ) 

DO 394 K = K I , KKF 
394 DX(K)=XH(K+1 ) — XH ( K ) 

AlPH1=AlPH*0.5 
39 DO 37 K=K I , KF 
37 B ( K ) = ALPH 

AP=1. J-( 1.0 /EXP (ALPH) ) 

A ( KF ) = AP 
Db 30 Mi-i=1,3p 
,< F 3 = 2 
,<ai = ka-i 

I F ( M‘-l- 2 ) 83 , 36, 3 7 
63 AlPH1=AlPH*0.5 
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n n 


ri= 1 

Ni'4 = 0 4 10=0 

Du 1 K = iK I » KKF 

DX ( K ) = X ri ( K+l ) — Xrl ( K ) 

1 H(K»M)=HGl+Cb*0.5*(XH(K)**2) 

H(KF»M)=H01 

AP=l.o-( 1 .u/EXP ( A lPH ) ) 

DO A K = K I » KF 

APC=AP* ( (HOl/h ( K.y.) ) * * 2 ) 

APO = 1 • O-APC 

5 P(<»M)=~( 1.0/AlPH)*ALOG( APO) 

WRITE(6»668) (K> P ( X , i'i ) , K = KI . KF ) 

DO 8 K= K I , KF 

8 P4(K)=P(K»M) 

HO(M) =Hul 
T ( M ) = 1 • u 
HE=H0 ( M ) 

I o=2 

CALL DVD 
CALL HDT 
PA I ( KF ) =0.0 
DO 6 K=K I , KKF 
J=KF-K 

6 PA I ( J ) = P A I ( J+1)-DX< J)*(DEN(J»M)+DEN( J+l.M) >*0.5 
DO A K= K I , KF 

4 VD( K ) = ( PA1 ( K )*B< K > ) / (H( K »M )**3*DtN ( K.M) ) 

CALL I N TEG1 (Xri( K I ) »XH(KF ) » 2 , VD , xF , VAlUE > IER ) 
EBF=VALuE-AlPH/ (.16 • 0* ( H ( K I »i*i)*PHZ ) **2 ) 

EF1=EBF 

DO 181 K=KA, KF 
DO 181 J = KA » KF 
BET ( K* J ) = o. 0 
181 CONTINUE 
GO TO 499 
56 M=2 4 GO TO 58 
5 7 Ki= 3 

58 NN=0 4 10=0 

IF ( MM . L E • 4 ) 5 8 J * 5 8 1 

580 Ho ( M ) =H0 ( Fi- 1 ) -20 • OE-6 
GO TO 983 

581 B i'i = M M 4 b H F = M F 

I F ( MM • L E . 1 5 ) 982.988 
982 RA1=RA* ( 1 .O-oM/BMF ) 

HO(M) = ( l.u-RAl ) *H0 ( M— 1 ) 

GO TO 963 

9 8 8 HO ( M ) = ( 1 . --C . 1 5 ) *H0 ( H- 1 ) 

98 3 DT(M)=H0(iV,-l J-HO(.'I) 

2 3 T(M)=T(i-i-i)+DT(M) 

I F ( MM . EO. 2 ) Go TO 161 

CALCULATION oF TrlE Initial OOtOoLD PRlGooRE oY linear 
EXTRAPOlAT ion. 

T 1 = ( T ( |V| ) - T ( 1-1-1 ) ) / ( T ( 0-2 ) - T ( i'i - 1 ) ) 

T2= < T < iV, )-T ( i'i— 2 ) ) / ( T ( 1 ) -T ( A~2 ) ) 

DO 751 K=KI » KF 

751 P ( K.M ) =7 1*P ( K.X-2 ) + T2*P ( K . to- 1 ) 

I o=2 6 11=2 4 CALL DVD $ CALL HDT 
14=3 4 CAlL INTEG 
HE 1 =H0 ( i v i— i ) 
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n n 


Hc=HO(M) 

GO TO 164 

161 DU 162 X = XI, XF 

162 P(K,M)=P(X,,v.-l ) 

XC = 4 6 w R I T E ( 6 * 2 0 u ) XC 
16=2.6 11=26 CALL DVD $ CAlL HDT 

164 EFl=EbF 

499 IF(NN+1.EQ.1 ) 1613. 118 

1613 DO 1614 X = X I , XXF 

1614 Pb(K)=PU»H) 

CALCULATION UF THE SYSTEM couATIuNS AND THE DEkIVATIVt.6 
OF THE SYSTEM EQUATIONS. 

118 I S=3 1 11=3 $ CALL DVD 
PNA = P(XA».'i) 

C6=EBF/(4.D*A(XF ) ) 

DO 71 X = XA * XXF 
Xx=X-XA+l 

CoDXP=C6/oX( is) * ( PI X+l ,M )-P ( X»H ) ) ■**■ D E)N ( n » i v l ) / V 1 6 ( K ) * { H ( X+ 1 « l»i ) 

X + h ( X » M ) ) ** 2 

E t- ( X K ) = C 8 D X P * ( H ( X + 1 » Pi ) + H ( X » M ) ) - ( P A I ( A + 1 ) + P A I ( X ) ) 

DO 71 J = XA * XXF 
JJ= J-XA+1 

IF(J.EQ.XA) 72, 78 

72 GoQ=G.D 

DO 73 I = X I * XA 

73 QUQ=QQQ+(U<X,I ) +U ( X + 1 » I ) )*P( I » M ) 

Q«Q=QQQ/P ( XA , M ) 

GO TO 74 

78 CmQ=Q(X,J)+u(x+I»J> 

74 R ( XX , JJ ) =-3.0*C6*CaoXP*OUu- (bc.T(X+l,J)+oET(K,J) ) 

IF(J.EQ.x) GO TO 1 o 

IFIJ.EQ.X+I) GO TO to 
GO TO '71 

75 SIGN=-1.U 
GO TO 77 

76 SIGN =1.0 

77 R ( XX , JJ )=R ( XX , J J ) + ( C6/ ( DX ( X ) *V I 6 ( X ) ) ) * ( ( P ( x+l ) -P ( X ) ) * ( D E N u 

X ( x »M) -DEN < X »M ) *VISD ( X ) /Vlo < x ) ) +6 I GN*DEN (mm) ) *.( ( H ( x+ 1 » M ) + 

X H ( X * M ) ) **3 ) 

71 CoMTINUL 

I E=XXF-XA+1 
467 IG=XF-XA 

C LIBRARY SubKOUTlNE, IN1SP* IS THE OP ERA T I OH uF MATkIX 

INVERSION .'.'HERE THE GAUSSIAN c.L I M I NAT I UN Ii'iElnUO IS USED. 

CALL IN1SP(R»IG»1.E-7,IEEk»51 * TEh.nI ) 

IF(IEER) 1503 , 15^3, 1304 

1304 WRITE (6, 200) IEER 
GO TO lou6 

1503 DO IGo XX=XI*IE 
A o = 0 • J 

DO 101 JJ=XI» IE 
101 Ao = AS+R ( XX , J j ) -*EE ( J J ) 

X=XX+X A- 1 
DP ( X ) =— AS 

P < K»M ) =P ( x»l-l) +DP ( X ) 

Po(K ) =P l x »M ) 

IvC CONTINUE 

/, k I T E ( 6 * 2 - J ) NN, Io 
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/■j K I T E ( 6 * 2 1 7 ) L F i 9 Lui* 

W R I T E ( 6 * 6 4 2 ) ( K 9 E E ( is. ) 9 K = R I > Ic ) 

\\A I TE ( 6 9 843 ) IKjDP(n) >,< = <A»^F ) 

Du 146 k = kA* KK.F 

— 1 • O UP ( Is, ) 

I F l SU ) 99^9 996» 146 

99 6 WK I TL“ ( 6 9 2 1 7 ) 6'J 6 GO TO 1000 
146 CuNTINJE 
13 3 P v. = 0 • 0 S Pu = w #0 
DO 106 k = is A 9 RKF 
P vi = P ,\'-f DP ( k ) 

1 ^ b P-j = PU+P ( K » i v '. ) 

P - O = A b S l P ;• / P u ) 

AAl=KA-l 

116 Ih(PRO-u.~jc2) lub » i w 89 117 

117 ,\.s = ,\N+l 3 IF (.\,\.LE«20> 111 9 192 

111 I h ( ,v i • E Q • 1 ) 1±29 1 1 j 

112 1 0 = 2 i 11 = 1 1 CALL DVj T CALL HjT 

PA I ( K F ) = u • u 

DO 791 .<=KI* KKF 
J = K. F - K 

791 PA I ( J ) = P A I (J + l ) -DX { J ) * ( DEi\ ( J 9 ; v t >+JL.\ t J+l 9 ■•■ ) ) *0.5 
CALL 1 i\ T L G 
Gw TO 1^6 

113 Io = 2^ 11 = 1 a CALL L V U 
981 C^LL HD I 

1^=3 4> ixO = 3 a CALL 1 1\ T 
GO TO 1^8 

158 IrtlO+l.E^.l) 1 1 o 9 15^ 

152 IF ( \>N+1 .EQ. 1 ) 1 53 9 1 1 8 

1 w 8 <j = 0 m \j i ^ 0 — v.. 1 . • , 

l\ 1 V = 0 

Du 2 24 k = is ,91 9 is A r 
p ;, = P + p ( K. 9 1*. ) -PS ( A ) 

224 Pv, = PO+P ( .< 9 .-. ) 

p r\ lor = A 13 S ( Pvv /P0 ) 

C Ti'L CALC Jl A 1 I Ji\ uh Oil I i\ T Lku^mIl o Imlli! Pkljojim:. 

i 3 2 2 fs — is I 9 A A 

3 2 2 V o ( K ) = P a I ( js ) / { n ( .x 9 ) * * 3 * J i_ ■ m ( is 9 i '- 1 ) ) 

0 .*n A = 1 . 1 •u/LX? ( alPH*P ( LA 9 .-, ) ) 

A ( K I ) = O . ’J <i> N. A 1 — s. a — 1 

D u 3 2 3 <\ — 1\ I 9 is A 1 

I 1 = K O I T I = is + 1 

CaLL I i \i T l 0 2 ( X n ( I T ) 9 X n ( I T I ) 9 1 9 7 U 9 rs \~ 9 V' A L O F. 9 I L K ) 

323 A ( I T 1 ) =t\ ( I 7 ) +VAL^t. 

Du 3 2 4 = » :s A 1 

U ! n — '..k' K A * a ( ,s ) / A { Is A ) 

•„* : \ = 1 • ~ — K 

I r ( -j K ) 3 3 1 9 JO 1 9 ^ d 4 

3 24 ? ( is 9 " ) = “ A l \J U ( 1 • 0 *.•« fx ) /ALPn 

G ^ 10 J J- o 

3 ji 1 is 1 = x ^ — 1 u /, i\ I T l. ( o ? ^ 1 / ) u ix 
l"A- 8 3c is — is I 9 .< A 1 

? 3 b P IF,'. ) = Pb ( ,s ) *P ( A A , . ) / P K A 

0 j o s I T h ( 6 9 3 4 4 ) ( n • rM is 9 1 -. ) 9 ,< = ?s I 9 is A ) 

?■ o Du 3 i , u <\ = .sl 9 is is F 

? .. 6 P c ( \ ) = P ( x 9 . J 

Ii ( '* " • 7 ^ • i *- ) 47 1 9 m-2/ 


60 



n n 


Ml IF ( PRO. i.E . 0 .00v,2 ) 197 » 227 

471 I F ( PRCJ# i_E • 0 • 001 ) 197* 227 

227 10=10+1 $ IF ( IO.LE.30) HI, 192 

153 EbF = (; . 5* ( CC ( KF ) -AlPH/ ( 16 . c* ( H ( K I »M ) *PHZ ) **2 ) +EF 1 ) 

EFl = CC(KF)-ALPH/(16.G*-(H(<I».'i)*PHZ)tt-*2) 

GO TO 118 
192 NU=1 

197 DO 193 K — X I * K.K.F 
SP = AbS(P(R»Pi) ) 

IF(SP-I.'j) 198, 19ti» 1000 

198 CONTINUE 

THIS CAuL Io FOR THE CALCULATION OF' LOAD bY THE dUdkuUTINE 
I N T E G 1 • 

CALL INTEG2(XH(i<I ) ,XH(KF ) ,2,P( 1,F(} ,KF, VALUE, IER) 
to ( M ) =2 . u* VALUE 
IF(MM.GE.7) 351, 632 

851 PUNCH 205, ( P ( K , M ) , K=KI» KF ) 

PUNCH 23 1, H ( ,< A , M, ) , TOM) 

WRITE (6, 668) ( K , PAI(K), K = KI, KKr) 

PRINT OoT-LuAD, CENTER FILM I il I CRN Luo , PRLbbUKL, 

852 A U T = A LP rt* w« 25/(H(KI» ,-i ) ) 

VIM) = — AP / ( C 4 * L 13 F ) 

WR I Tb ( 6 , 1 - j 9 ) AuT 
V C = A P # ( HQ ( M ) **2 ) *PhO/ALPH 
W R I T E ( 6 » 1 u 0 9 ) VC 
DO 853 K = K I » KF 

853 VO ( K ) =XH ( K ) /H ( K *M ) **3 
DP ( K I ) - j t-j 

DO 854 K=KI,KKF 

854 DP(K+1)=DP(K)+0.3*(VD(K+1)+Vd(K) ) #DX ( K ) 

VF = — AP / I C4* I DP ( K ) * AlPH— ALPn / ( 1 b • O'* I n ( A I , ) *Prl2 ) ** 2 ) ) ) 

Wr< I Tt ( 6 » 2 1 7 ) VF » DP ( ieF ) 

WK I TF. (6,220 ) .Mi'i , << ( n ) , no (,>•,) , f ( i-i ) , D T I-.-i ) , V ( n ) 


WRITE (6,210) 

( K 9 

P ( 9 1M ) * 

i — i 

,/ 

it 

/ 

9 fxF ) 

WRITE (6,211 ) 

( K 9 

h ( k, 9 /i ) ^ 

K - \ I 9 

kF } 

WRITE (6,215) 

( K 9 

D { rC • : v i ) 9 

V 

II 

-<s F ) 


I F ( Mf- • L L • 2 ) GO To 86 3 
DU 861 K = K I , kF 

8 61 VD ( K ) = V ( r-i ) * ( oi-lEG 1 ( i-i ) *0 ( ,< , r-i ) +Oi-iCij 2 ( i v l ) *D ( il , ,-i- J. ) + UH Lb 3 ( H ) *D ( iv *r'i-2 ) 

X - 1 . • ) 

W R I T E ( 6 , 2 0 9 ) ( K , V D ( K ) , K = K I , KF ) 

8 6 3 WS=4 • u# ( ( PH/t ) ^ ^2 ) ( i'i ) 

to T = 1 • 5 * w ( i-i ) 

WRITE ( 6,21 7) T , ,mS 
Iu=2 i 11=26 CALL UVD $ Cm l L bDT 
557 IFIMM-2) 50, 30, 32 
52 Du 41 M =1,2 
T(M)=T(H+1) 

HU ( M) =HU(|V;+1 ) 

DO 41 JC = K I » KF 
DEN ( K , i’i ) = DE i\ ( K , Pi + 1 ) 

D l K , M ) = D I K , ;«i+ 1 ) 

H ( ,K » M ) = H ( rx , Pi + 1 ) 

P ( X * X 1 = P ( n , ,v, + l ) 

41 CuNT INOE 

IF (HOI M) . lL . lC.uc-6 ) 1 ;;cu » 1099 
1099 IF (MM.EU.3 ) 317, 316 

317 ‘M 3 - ,v ( i-l ) i H 3=HO ( i- 1 , ) i T3 = T(i-i) + V3 = V(Ri) 4> dT3 = dT(.-i) 6 uu Tu 30 
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3 1 a '/j ( MM ) = o ( M ) ^ T ( , v ii‘i ) = T ( .*i J ^ rio ( : -u v i ) =ho ( i-i ) ^ V ( .-ii*. ) = V ( i-i ) ^ o T ( i*u v i ) = o T ( r*i ) 

^ CwNT 1 N’Jl 
i w'u I f t. ( 6 * 6 v J 
1 j ^ J - V * 

U j 5 8 3 I - K 1 » k.'i 
I F ( I • E 0 • 3 ) 519 9 5o3 

519 W( I )‘ = >v3J>H( I )=H33 7 ( I )=T3$V( I )=V3£JT< I )=DT3 
583 WkITE(6*221> I, * ( I ) ,HQ ( I ) > T l I ) , 0 1 ( 1 ) * V l 1 ) 

STOP 1 

2-;0 F0RV.AT(M5® 2F10.3, FlC.l* FiO.O) 

Pul FORMAT < 4F 1 5 • 8 } 

2 w 2 format ( 5 E 15.2 j 

2 v 3 F o R M A T ( F i 5 • 2 ? F I 0 • 2 , 515*2? r 1 u • 3 ) 

2^4 FoRMAT (2Flu.l, lKi^j F 1 G ■ i ) 

2 G 5 Format ( 6E12. 5) 

2w6 Format ( 6EI2. 5 } 

2^7 format ( ofiu. i j 

2 v 9 FORMAT ( jX//hCX» * RELAX I V !_ \/ E L oC I f Y , V ( X 9 . v i ) = * / / 2 X , * is * » 2 j X » * n* ,2 G A 9 

X * K* 9 2 o X 9 * k* 9 2 X 9 * X * 9 2 0 X , * X * / / ( b ( I 5 , L 1 o • 7 9 J> A ) ) ' ) 


21 c 

format ( 3x//b ox , 

* P ( M 9 M ) 

= * / / 2 X 9 ■* 

N* 9 

20X , 


*N 

* 9 

20X9 


X2 - X » *K*, 20X » 

- k 'M2 0X i 

* <“*//( 6 ( I 3 9 

L 1 5 • 

7 

, 

3X ) 

) ) 

211 

format ( ax //box , 

* ri ( k 9 i-i ) 

=*//2X9 * 

,N* 9 

2 0X9 


IN. 

* 9 

2 0 X 9 •* k * 9 


X2-X? K 9 2 0 X 9 

^ i\'M2 0 X 9 

■*R*// ( O ( 

I 3 9 

Li j • 

/ 

9 

^X ) 

) ) 

21b 

format ( jx//^ jX > 

L/ ( .^ 911 ) 

= * / / 2 X 9 -* 

N* 9 

2 0X9 



* 9 

2 0X9 -*K*9 


X 2 - X 9 * rC * 9 2 0 X 9 

* IN * 9 20X9 

/ / ( O ( 

13, 

Li J • 

/ 

9 

^X ) 

) ) 


212 F 0 RM 1 AT ( o X / / 1;A> * ] n t r\ I G i u o v o Y R K l o o u R l 9 F ( i\ 9 .- 1 ) = * / / 2 X 9 



X*K* - 

9 2^ 

X 9 

* 

In * 9 

2 ^ A 

9 

* 

A'* 9 


2 

ox 

9 

'X- lN -x- 9 

2 

O A 

9 -* * 

» <i 

0 x 9 

*rx*// ( 0 ( I J » 


XE 1 3 • 

.79 

) 

) 

) 
















22 

3 FoKr 

'’AT ( 

Ini 

9 

3X 9 


= * 

9 

13 

9 


^ A 

9 

-*V\ ( i'l 

) 

— - j'. 

9 1 1 5 

. /, 

X 

*v»o=* * 


XB 1 6 « 

.79 

sX 9 


■*riO { 

Mi ) = 

'* , 


£13 

• 

7 

9 

■jX 

9 -* V 

= 

* 9 

11 

A 
~H 
1 J 

7 > 



22 

0 F 0 R ( 

•AT ( 

1H1 

9 

3 X 9 


= 

9 

I 3 

9 


5X 

9 

* ( rT 

) 

= -x- 

» Elb 

. 7 , 

3 X 9 

*H J ( i». ) -* , 


XE 1 5 < 

.7, 

3X 9 


* T ( M 

)=* 

9 

L 

15. 

7 

9 

3 

X 9 

"* 10 T 

( 


= * > b 

. 

i.O 
■ — t 

7 9 2 X 1 

> *V(ii)=*»Elb.7) 

22 

1 FoRf- 

'AT ( 

3X 9 

-K- 

1 *'i — ** 9 

13, 

^ A 

9 

*w ( 

, ‘i 

) 

— -X- 

9 E 

13.7 

9 

3 X 

, * ri w ( 

.'i ) = 

* 9 



X E 1 5 • 7 9 3 X 9 * T ( i-t ) = * 9 L 1 5 • ( 9 3 X » * ^ T ( i-i ) = * 9 el 1 3 • 7 » 2 X 9 * V ( 1 ) = * 9 b 1 5 • ( ) 

6 vo FoRKA f ( jX / / / j I 3 9 GX 9 -*EP Io l oo oMALL*) 

231 FORMAT (2L2o.io) 

6 A w F o K M A F ( 3 X 9 ■* 1 * — i < 9 PklooJRL o A I A ™ / / / ( o ( 1 j 9 F i o • b ) ) / / / ) 

7 G ^ format (^ x» !:✓//(<+{ i^, fia.oh) 

7 1 J F ^ k ■■-! A T ( 5 X 9 I o 9 A r 12*6} 

7 31 F O k X AT ( ini » 3 a, *Fi = * > 13 ? 'J a 9 * « ( . ■ - J = # 9 l! j» / j 3 a 9 *>ro = * 9 L 1 3 • / 9 

X 3 A 9 *HO(-i>=*» L 1 3 • 7 ) 

666 format ( els* 7Ei6«y; 

669 FORMAT ( 1H1// (6( 14, IX* E16.8))) 

668 FoRMAT < 2X// < 6 < 14 9 IX* FI 6 . 0 ))) 

6-3 FoR^AT ( 3 X// 5 JX , * 5 A 1 R i X r(.^»o)=*//(6(2I4,E12*:>) ) ) 

1 %. J d F o- K MA T ( 1 H 1 / 2 X 9 -•*'■ V = * 9 L 1 6 • d ) 

1 C 9 F v k A T ( X 9 -* v 1 = * * ii 1 6 • 6 ) 

16 - C F oi< A f ( 4 5 X 9 *,i(M)=** El 6 . 3 ) 

l.io Format (4:>x, *trjc v = *» eio*6) 

2 1 o F ^ aMA I ( o ( F i 1 • 0 9 F 1*9 o ) ) 

217 FoR.iATt /) 

3 2 9 F ^ i’i ■ ; A 1 ( <_ X ? * 1 ■ . A T k I X 1 j SI Ko ^ L n i\ ? 2 X ) 

0 A 4 FvK.-'Al ( * 4 o X 9 ^ P ( ;< 9 : ■ 1 ) = »" / / (61 I 4 t 9 2 X 9 c i o • ( ) ) ) 

8 41 FvK/.A T ( ^ ^ a 9 «■ ,\ i 1 ( n ) = * / / ( 6 ( 1 h»^a>h^*j) ) ) 

84 2 F o R M A T l -r o X 9 -»- __( \ ) = «//( 6 ( I 4 ? ^ X 9 1 : ^ ^ • 7 ) ) ) 

8 F" w R AT ( X » uP ( k )=■'<■// ( 6 ( I 4 9 2 a 9 E i 3 • 7 ) ) ) 
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non n, n 


I 


C^LCuLAT I Jr* ut F I L'*i Tn ICR NE Go A : -4 l> CucF r 1 C I c.i\ I G Or LAG RANG I AN 
T 1 iKLL — Pu I 1 \ T uoADr<A 1 uR«_. 

C ^ r. V.U N / A / x F 9 A l P n 9 t_ N , E 0 9 i 1 E 9 m E 1 9 A L P H 1 9 . w 'u v i 

COi-ii-:UU I 9 KM J N i\ 1 > KKp ) is. r 2 9 i\H J J V* 9 nn c. 9 I O 9 j >i 9 N K 5 Po * N * . v * 9 

1 P 1 • 1 9 PC* no* uP * E o » >« 9 v'v >» * 1*1 * Km * kU * P in k » P U i 9 X ri * OX * w I 9 12 * 

2riG*H*P*Vlo*v’lou>9L^i_r'j*UiiND*T 9 o 7 * V u * uP * H \J * o l. * a 9 0 * C * G u * m«*'i * P G * 

3Ci‘iEG3 9 u« v ic.G2 joi'IEnjI » [)u> » l> j * C 4 * C 2 * i I > ^ » Pa i »clT 9 P 1 * P 2 * v 1 G 1 * 

4 V i o 2 * V I 3 *0PV*C6 *LdF 90 LJBV 9 C C * L c. 9 P n a * D P R * R 
u 1 . v lEi\iG 1 oi\i Xri ( 2 1 ) * DX ( 2 1 ) 9 0 t 2 i * 3 ) 9 n G ( 2 i 9 3 ) * ri t 2 i * 2 ) * u I (jJ.)* 

I P ( 5 1 * 3 ) * V 1 15 ( 3 1 ) * >/ 1 G u l 3 1 ) »BEN(51»3) 9 12 l N D ( 3 1 9 2 ) * T ( 3 3 ) * V ( 2 3 ) * 

2DT ( 35 ) * VDl 31 ) *L>P ( 51 ) *nQ( 30 ) * » ( 3 3 ) *0 1 ( r 3 i ) * A ( 3 1 ) , b ( 3 I‘) * CC { 3 1 ) , 

3F l 3 1 ) 9 R ( 3 i » 2 i ) >ocT ( ji *31 ) * ju l 3 i ) *Po(3i) * wniiG2 C 3 ) > G'i*iLG2 ( 3 ) * 

40KEG1 ( 3 ) .DENT (31) *0(31*51) * P A I (31) *oubV ( 31 ) *C(2*3*3l) *EL(31) 

DO 30 J .< = n 1 » i<F 
D 3 — \j • o 

DO 3C1 J=n I , <F 
3w 1 Do = Q ( R * J ) *P ( J *M j + D6 
3 u 0 D ( K * M ) = - C 6 * U G 

DO 303 ,n = ,<I* <F 
I F ( KF3 • EQ • 1 ) 22* 42 
42 IF(M,V|.E0. i) 21, 2 2 

21 Fib ( K * iA ) = rlG ( n ) + J • 3 * C 3 * ( X H ( »\ ) * * 2 ) 

Fi ( !C * i v . ) = i~j G ( K * 1 ; i ) + O ( K*n) 

Gw TO 303 

2 2 H G ( K * l v l ) = H 0 ( . • 1 )+ '.j • o * Co* ( X H ( n ) * * 2 ) 

Fl(K* Vt ) = HG(K9 ) +D ( fs 9 r'l ) 

3 w 3 CONTINUE 

HO ( M ) = H ( KF * F; ) 

341 I F ( M M — 1 ) 2 IO 9 3 1 U * 2 il 
311 DT (Ri)=HJ(n~l )-H0(n) 

T ( y. ) = T ( .*;— 1 ) + D T (,M) 

+ F ( — 2) 31 w * 3 2 2 * 3 3 1 

333 0imEG1CM) = 1«^/DT(.v,) 

Gw TO 310 

331 Cv.EG3 l.-i ) = ( 1 ( r.)-T ( ,-.-1 ) ) / ( ( T ( ,-,-2 ) - 1 ( »* 1 ) ) * ( T ( ,-.-2 ) -1 ( n ) ) ) 

Qi*.EG2 ( 1*1 ) = ( T (.-.)- T ( 2 ) ) / ( ( T ( 1 )-f (r-2) ) * ( T ( ,-i- 1 ) - T (.'»))) 

On EG 1 ( M ) = ( ( T ( .'i ) - f ( n-1 ) ) + ( T ( n ) - T ( r.~2 )))/(( T ( n )- T ( ■'*.- 1 ) ) * ( T ( n )-T ( .-.-2 ) 

X) ) 

310 CONTINUE 
RETURN 
E,A) 

subroutine integ 

CaLCuLA l I u.\ w r Inc. I iU EG^al A 1 4 u UlKiVATiVl w f ■ l ti l- i 1 ^ \ u K /-\ L 
1 1 m H I C ri T n u 1NT coKanj I G T I i^lRIVmT IVl uF T n E 
P R 0 L) U C 1 0 F 0 ci\3 i T Y /\hLi F 1 L n l n 1 LkNLoo# 

C o f «1 |V 0 1\ / A / i\F * AlPH * iiiN? Eu* liii* mliI* mLPI'iI* . v h v 1 

Cwi v l f 'OK in 1 * nA * nnI 9 nnF , i\r ^ 9 n P d 9 V 9 n i\ ii 9 io » \ 9 i\l\ 9 PO 9 n 9 n> 

1 P i v - 9 P C 9 N 9 Lp 9 L • J 9 A 9 A 1 A 9 /* .‘I 9 i\ A 9 KJ 9 P N K 9 P 0 I 9 X M 9 U X 9 O 1 9 9 

2 H O * H 9 P 9 V I J S V I J J 9 U C . ^ • D L iNJ D 9 I 9 \J T 9 V 0 9 J P 9 r! O 9 O L 9 A 9 D 9 C 9 O 0‘ 9 A i' ' l 9 P G 9 

30UEG3 * 0:*il.o 2 * O t*i ii G 1 *^i_i\T j'jjjJojCAjO > I I * u * P / \ I *uiil *P1 *P2*V I o 1 9 
4 V I 32 * V 1 1 ? 3 * D P V 9 C 6 9 ibF * o J i j V 9 C C * L E 9 P K m 9 D P R 9 R 
D I ENG I OJs! X hi { d 1 ) 9 D A ( 3 1 ) *0(31*3) 9 ri G ( ’J 1 * 3 ) 9 ri ( J i 9 _> ) *-u I ( 3 j. ) 9 

IP ( 3 1 9 3 ) 9 V I G ( ul ) 9 VI ( Jl ) , Dl-im ( 3 1 9 3 ) * JEN0 ( 3 i * j> ) » T ( 3 3 ) * V ( 3 3 ) * 

2D [ (33) * V 0(31 ) 9 DP ( 3 1 ) 911 OOG) * w ( :> 3 ) * o L ( 3 1 ) * A l 9 1 ) * b l 2 1 ) * Co ( 3 1 ) * 

J F ( 3 1 ) * R ( b 1 * 2 1 ) 9 OC.F ( 2 I *31 ) * w i3 ( 2 1 ) * P w ( 2 1 ) * w.*i n ii 3 ( 2 ) * oriuvj2 ( 3 ) 9 

4J, v i i: C j 1 ( 3 ) * -j c . m i ( 2 I) 9 xji ( 2 1 * 2 1 ) 9PmI(2x) 9 o ^ w V’ ( 2 1 ) »U 1 j .3 » ) j 11 ( j i ) 

< A 1 = K A - 1 

22 I F ) 3*1* 2 
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1 i) w b .< = i\I ? k r 

D L i\ T ( f\ ) - ( i; ( \ j i*i ) " Jl is ( i\ > .'i~ 1 ) ) / U T ( r. ) 

sJL.' I K. ) — i \J ( .'s. 9 i '| ) i_J ( i'n r-\ 9 i‘ I ) ” 1> ( in 9 ■'! 1 ) l in A 9 i l . “* 1 ) ) / L ) I ( i*l ) 

b O L. ( K ) = [>’ nt’i l *n 9 tO ( i\ ) — 1 • \J ) + i 1 ( .n. 9 .*i ) “ r 0 C is 1 ( in ) 

P A I ( K F ) • j 

D 0 4 N = .‘v I 9 In F 

J = KF-;< 

4 PAI ( J ) = PAI ( J + l ) + i';X ( J ) « ( SE ( J ) + SE ( J+l ) ) #0. b 
Gw TO 31 

2 jJ >J 1 \~K i > N F 

0 'iL M T ( A ) - Ji-ltl'C J ( -I ) ( A 9 ■ -I — 2 ) + , -j L (j 2 ( 1-1 ) '"'JLi'J ( A * r'l — 1 ) + 

X 0 . ■ . c. G 3 ( « l,- ^ ( , v 9 - i ) 

5 ^ ( K ) = J.-l L 0 3 ( .1 ) & JLis ( ,n 9 1*1 — 2 ) * ' 0 ( A 9 *•<■“ 4 ) + «J:-I iL sj 2 ( -I ) Jc.i'j ( fN > .'i— 1 ) * J ( K 9 fi — 1 ) 

X 4- • i b C j 1 { i * i ) '« r. i'i ( >\ 9 it ) ^ \j ( is 9 ) 

10 5 u ( K ) = “ J L l\ ( K ) I' I ) + O i> ( I\ ) '*" H G ( K 9 1*1 ) L* L. I * 1 ( fN ) 

P r-\ I ( K F ) = J • ^ 

DO 2 b ,< = .vl ? kkF 
J = K F' - K 

2 b PA I ( J ) =PAI ( J + l ) +uX( J )* ( S£ ( J ) + SL ( J + l ) )*0.b 
31 IF ( IS.Ew.3) lb? 20 
lo D-j 16 K = i< A 1 9 KF 
Li Li T ( K F 9 a ) - o • -j 
16 CONTINUE 
in A 1 — K, A 1 
K F 2 = K. F - 2 

0 ^ lb sj K = In A 9 In F 

[)w 1 D J J = i\ M 9 t\F 

|\ ( In > J ) - w * . 4- IF ( J • C U • \ A ) i J 1 9 ±32 

1 b 1 1 b 3 i=;sl > n A 

1 b 3 R ( K 9 J ) = K { In 9 J ) +0 ( |'N 9 i ) * P ( I 9 .■'! ) / p ( fN r\ 9 > v i ) 

G T 0 1 0 j 

1 b 2 R ( K 9 J 5 = vrf ( 4 9 J ) 

15 0 COM T I Fi J 2 

Do 1 0 1 tn = n A 9 k F 

Dj- 1 f D 1 J = .\A> ,0F 

6 0 = w • - 

D'w 1 0 r-j = io in i\ F 

1^5 56.b=S65-U2N ( r-j 9rl ) 9 J ) * ( Xu ( u+1 ) - Xu ( u-l ) ) 

6 N -' 1 = S s S + l) E i’-J ( 9 ' v • ) ' x * ( Tn 9 J ) “rt'f)X ( 1 ) 

5 o = S S - 0 n i vj ( <\F 9 .-I ) * rt ( A F 9 J ) X ( in in F ) 

hi L T ( K 9 J ) - 0 . v i L.b 1 ( . v i ) 0 6 ^ w 5 * 0*5 

1 F { J - ,< ) 1 vJ 1 9 1 1 ; 7 9 1^6 

1 J 7 bLT(K»J)=i^LT { X • J ) + J c. ^ i.> ( J 9 i v i ) OX ( J ) * ( 1 6 1 £ .'*■ ) ( J 9 ) — 1 • Q ) * 0 • b 

GO TO 101 

1 .. o : j z_ T ( K 9 J ) = ui ii T { is 9 J ) + iJ E .m 0 ( J 9 r*i ) ~ n ' ( X m ( ^ ■+■ 1 ) ~X hi ( 0 — i ) ) ^ ( Ori c o 1 { .*i ) 

X -* n { J 9 j *1 ) - 1 * ) -* 0‘ • 5 

1^1 C-MT 

D ^ 3 3 < = \ A 1 • .< F 
3j> j T ( < F > i\ ) = v • 

G w T 0 6 ^ 

j D 4 ; ; X = A A ? fN F 

U ^ ■. J — i'N 9 i\ F 

j ■_ T ( K ? J ) = • s,' 

4j c^n\ju 
o 0 D ^ 4 / in = -N I ? ‘n F 

4 7 V- : ( K ) = D C A } ^ PA I ( x ) / H ( k , ’ ) » 3 ;; 'D Fi - « ( k ? ) ) 

Call i ,*s r £ g 2 ( x r i { .\ i ) 9 x n t : s r ) 9^9 v o 9 nF 9 v alj £ ? 1 £ x ) 

U £ nF i — V L_ vy iz 
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n n 


! 

2 0 C^H 1 I NU c. 

RETURN 
E NO 

b 03 R 0 J T I i\E 0 V L) 

C.^LCuLATiOiM «jh J l U ^ i T Y j VIoCuoi i Y a l\ 0 1 1 1 c. 1 r* lJlK I VhT I Vlo 

WITH REoPCCT TO PREoauR E. 

C wMHOM / A / ;\ F * A l P f i 9 L <\ 9 ED 9 f 1 1_ 9 i i 1 9 aLPhl » ri. v i 

Co 1V0.M n I 5 \A» \\I > R k F 9 Kr 2 9 kF d 9 v » KK.E 9 1 o > *"m 9 Ni\ 9 PO 9 R. 9 n 9 

1 P. v i 9 PC 9 Is o 9 t.P 9 il J ? 9 '.V 9 i v i 9 r;A 9 i\U 9 PhK 9 PU ] 9 Xil » JX j D I 9 u 9 

2 H Lj 9 H 9 P 9 V I 0 * V I o 0 9 U l_ N 9 0 c. N !.) 9 1 9 J I 9 V l> 9 Ly P 9 ri CJ 9 O Q 9 A 9 LJ 9 C 9 0 u 9 A 1*1 9 P b 9 
.sCriEG 3 9 u**‘E.G 2 9 j.mlu 1 9 u li\ 1 > Ou sCA j 9 i I » j 9 P A I 9 o L T » P 1 9 p 29 ‘VIbi 9 
4 V 1 5 2 9 V I . j 3 9 Lj P V 9 C 6 9 E to F ? b ! J b V 9 C C 9 L l 9 H K A 9 0 P i < 9 R 

D 1 M E Mol u N X ri ( 3 1 } 9 J a ( 3 1 ) 9 D ( 2 i 9 j ) 9 n 0 ( b i 9 3 ) 9 rl ( o jl 9 ) 9 O I { j i ) 9 

1 P { 3 1 9 3 ) 9 V I O ( 3 i ) 9 V I o ( 13 1 ) 9ijLli'4(Di93) 9L^l_l\0(3>i9J) 9 T ( Jj 3 ) 9 V t j> v j ) 9 

20 ^ ( 3 u ) 9 v‘ u ( 3 1 ) 9 iO P C l? I ) 9 1 1 G ( I? u ) 9 vj i i 'j ) 90 t ( 5 i J 9 A ( 'j 1 ) 9 u ( b 1 ) 9LL ( 3 1 ) 9 

3 F ( 5 1 ) 9 R ( b I 9 b 1 ) 9 l:»E1 (51*51 ) 9 bL> ( 5 1 ) 9 P 5 ( 5 1 ) 9 ui v i L G 3 ( 3 ) *9 0*1 c b 2 ( 3 ) 9 

4 0 1 *. L G 1 ( 3 ) 9 0 E 1 m I ( b i ) 9 w ( 3 1 9 5 1 ) 9 P A 1 ( b 1 ) 9^ubV ( 3l ) jC ( z » j »3l ) 9 L E ( 5 i ) 

Ir(IS-2) 9 2 o 2 9 2 0 7 

2 - C VI S < K ) = LXP ( A L P H * P ( 9 ■*, ) ) 

V I SD ( < ) = A l P H #V I o ( .<) 

DlN ( K 9 /. ) = i • 0+( E.\»P (,<. , ,/))/( i # c+EL)*P ( R 9 H ) ) 

DLND( k 9 . 1 ) = E K / { ( 1 . w + P i A ) *Lu ) **2 ) 

Go TO Zo 

2 l J 2 0w 2 g3 n> = KI 9 KF 

VI 5 ( K ) = l A P ( ALPHtfP ( k 9.v ) ) 

V I bl) ( K ) = AlPh*V I o ( k ) 

0 L K ( K 9 >| ) = I • w + ( C. I'M * P ( N 9 I-I ) ) / ( 1 • + L U * A ( iN 9 1-1 ) } 

2-3 Dr. NO ( X 9 .-■, ) = L A / ( ( l . -y + P ( K 9 0 ) * c. L ) * * 2 ) 

GO TO 2^6 

2 ^ 7 DO 2 0 V i\ = fN t-\ 9 K .< F 

P3= ( P { K + l 9.-: ) +P ( K 9 r, ) ) / 2 • 0 
VI 0 ( K ) = E X P ( ALPn*P3 ) 

VI 00 ( •< ) = J . 5*ALPH*V I S ( <) 

J i F ( K 9 >i ) = I • ’w + ( L*_ u « P 3 ) / ( 1 • 0 + l u * P ) 

2 - V 0 l ; N I 0 ( .< 9 ••■ ) = j • 3 * E i\ / ( ( 1 • + P 3 * E J ) ■* * 2 ) 

If (II*"w*i) 46 ? 2 0 

46 ,< A 1 = K A - i 

1 F ( f \:A + 1 • ..w • 1 ) 2-9 03 

Ov 4 7 i< = ^ I 9 kA 1 

47 p ( Ks 9“-' ) =P ( .N J ( N A 9 . ■ ■ ) / D < A 

2 - CjM TIMJl. 

Rl- TUPi\ 


?c 
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C WEIGHTING FUNCTION PASSED TO THL .-lA I N PROGRAM 

FOR STORAGE. 

I ERR RESULTANT ERROR PARAMETER. 

REQUIRED SUBPROGRAMS - NONE 
COMMON STORAGE - 

TnE WEIGHTING FUNCTION C IS STuRcD In THL HA IN PRoGRA.*! AND 
R L 0 U I R L G 1 H L FOLLOw I Nb DIMENSION 3 T A 1 1 M L N T ft H c R E u . G L • N P • 

ERROR INDICATIONS - 
I ERR = 0 INDICATES NO ERROR. 

ILRR = 1 INDICATES NP IS lEgS -TriAiM A. 

I ERR = Z INDICATES THE LIMITS OF INTEGRATION AkE uUT OF 
THE RANGE OF THL TABlE. 

COMMON/ A/^F-t AlPH» EN? EL' > Ho HEl* ALPnl, MM 

COMMON A I > K.A » KKI? <RF * i<F2? KF3 > V» K K £ * IS? N? NN* PG* & ? M? 
1 P i v i ? P C 9 K S 9 E P 9 nU » ft ? '.v v-J 9 •■>'* M 9 R A 5 R 0 * P i\ R ? P D I * X h * D X * u I 9 0 9 

2Hb?H?P>Vlj>VIjDjLjEi\jDLND>T ? D I )VL , »L)P)HiJ>ot:>AjD>C>Si)>A »'i 9 P S * 

30MEG3 9 OMcG 2 ? OMLG i » OLNT }Du>Uo jC4»Cd » I I >w>PAl »olT * P 1 ?P2*VIS1* 

AVI 02 ? V I o 3 »UPV > C6 9 ELF ? u U 3 V »lC ? E c ? P K A ? D P R ? k 
DIMENSION XH ( bl ) 9 DX (31) 9 0 { d 1 »3 ) »HG( bl »3 ) »H ( 31 5 3 ) 9 Jl ( 31 ) > ~ 
lP(31»3)tVlS{al)»VIoU(Dl)»0c:N(Dl»3)»uLUU(3l»3)tT(33)»Vl33)» 

2DT ( 3b ) »VD( 31 ) *DP ( 31 ) *H0 ( 30 ) »/;( 35 ) »SL ( bl ) *A( bl ) »b( 31 J »CC( 31 J • 

3F ( 3 1 ) » R ( 5 1 9 5 1 ) 9 B E T (31951 ) ? oD ( b i ) ? P S ( 3 1 ) 9 UM E G 0 ( 3 ) 9Ui v .LG2(3) » 
40MEG1 ( 3 ) 9 OENT ( 31 ) 9Q( 31 9 5 1 ) 9 PAI ( 51 ) 9oJbV < 31 ) 9 C ( 2 9 3 9 51 ) fEE ( 51 ) 

N P MUST BE GREATER THAN 3 
IF (NP.lE. 3 ) GO TO 96 
CALCULATION OF INTERVALS OF X 
Nh=NP-l 
DO 10 I =1 9 N H 
10 DX ( I ) =XH( 1 + 1 )-Xh( I ) 

DO 2 0 1=1? N H 

DcFINt uOiiFh ICI £NTo OF FIRgT PaKAdOLA 
IF ( I . EQ . 1 ) Gu TO 13 

C( 1 9 1 9 I ) =-DX ( I )**3/ ( 6 .0*DX ( I -1 ) * < L)X ( 1-1 ) +DX ( 1 ) ) ) 

C ( 1 9 2 9 I )=DX( I ) *( 3. j*DX ( I -1 )+Da ( I ) ) / ( 6.0*QX ( 1-1 ) ) 

C ( 1 9 3 9 I ) = DX < I ) * ( 3 .0*DX ( I -1 ) +2.U*oX ( i ) ) / ( 6 .0* ( UX { 1-1 )+uX ( I ) ) ) 
lb CONTIN E 

IF ( I • E Q • N H ) GO TO ^0 

DEFINE COEFFICIENTS OF SECOND PARIdOLA 

C ( 2 ? 1 ? I ) =DX ( I ) * ( 2.0*L/X ( I ) + 3. u*DX ( 1 + 1 ) ) / ( 6.u*( uX ( I ) +DX ( I + 1 ) ) ) 

C ( 2?2 9 I ) =DX ( I ) * ( DX ( I ) + 3 . 0*oX ( 1 + 1 ) ) / ( 6.0*oX ( 1 + 1 ) ) 

C ( 2 » 3 9 I ) =-DX ( I )**3/ ( 6. C*DX ( 1 + 1 ) *( DX ( I )+DX ( 1 + 1 ) ) ) 

20 CONTINUE 

EnTRY I N T L G 2 

I N IT I AL 1 Z E S^MMA I I On VARIABLE 
VALUE — 0 . u 

I r ( G 2 — G I ) 4G?92?3u 
B IS G R 11 A T L R I H A N A 
30 AL I M=G1 
B L I M = G 2 
SIGN = l.j 
GO TO bu 

A IS GREATER THAN B 
40 A L I Y = G 2 
BL IM=G1 
sign =- 1 .: 

bO N n = i\ P — 1 

IMKCT.EG.1) 12b* 123 

Calculation of integral uvlk gob interval 
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123 DO 80 1=1* Nh 

5 u b v ( i ) = o • j 

IF (XH( I ) .EQ.ALI.-1) SubV ( I ) =C ( 2 . 1 , I ) *F ( I ) +C ( 2 » 2 » I ) *F l I +1 ) 

X+C. ( 2 » 3 » I )*F< 1+2) 

IP < XH( 1+1 ) .EQ.bLlM) bU3V < I ) =C ( 1 , 1 , 1 )*F ( I -1 )+C ( 1 ,2 , I )*F ( I ) 

X+C ( 1 *3, I ) *F < 1 + 1 ) 

IF ( XH( I ) .GT.ALIM. ANJ.XH ( 1 + 1 ) .LT .bLIrt ) SuBV ( I ) =0 . 5* ( C ( i * 1 » I ) 

X*P (I— 1)+(C<1»2»I)+C(2»1,I ) ) # F ( I ) + (C(l»3*II+C(2»2*I) J * F ( I + l) + 

XC t 2 » 3 * I )*F( 1+2) ) 

80 VALUE =VALUE+3UBV ( I ) 

VALUE=5IGN*VALUE 
go TO 92 

125 DO 110 1 = 1* 'NH 
SUBV ( I ) =0.0 

IF (XH( I ) .EO.ALIM) 111* 110 
111 IF(I.EQ.Nri) i 1 3 * 1 1 A 

113 SUBV ( I ) =C< 1 * 1 , I )*F l'i-1 )+C ( 1 .2* I) #F ( 1 )+C( 1 ,3 , I )*F ( 1 + 1 ) 

GO TO 120 

114 IF(I.GE.2) 115, 116 

115 SUBV ( I > = 0 • 5 * ( C ( 1 » 1 * I > *F ( I -1 ) + ( C ( 1 » 2 » I ) +C ( 2 » 1 » I ) ) #F ( I ) + ( C ( 1 » 3 » I ) 
X+C ( 2 *2 , I ) ) *F ( 1 + 1 )+C ( 2 ,3 , I ) *F ( I +2 ) ) 

GU TO 120 

116 5uRV( I ) =C ( 2 » 1 , I ) *P ( I ) + C ( 2 , 2 » 1 ) *F ( 1 + 1)+C(2,3»1 ) *F l 1+2) 

120 VALUE=S !GN*SUbV ( I ) 6 GO TO 92 

110 CONTINUE 

* SET ERRUR PARA Me. TER FUR ,\iORMAL RETURN 
9 2 I ERR = 0 

RETURN 

* SET ERROR PARAMETER FOR TOu FE ^ POINTS 

96 I ERR = 1 
RETURN 

* SET ERROR PARAMETER FOR A AND/OR b OUT OF RANGE OF TAoLE 

97 I ERR = 2 
RETURN 
END 

* END OF RECORD 
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LIST OF SYMBOLS 




Half of Hertzian width 
coefficient of density 

coefficient of density 

constant in deformation formula 

constant in deformation formula of cylinder 1 

constant in deformation formula of cylinder 2 

coefficient of deformation formula 

Deformation 


Equivalent Young’s modulus 
Young’s modulus of cylinder 1 

Young’s modulus of cylinder 2 

Film thickness 


Rigid center film thickness 
center film thickness 


geometrical film thickness 
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V -t 


w = — 
m R 


v-vv 

1 < -V Z J ) 


P=I 


p = — 
HZ E 


“k.J 


v 

T -r t 


Minimum film thickness 

A dummy index 
See Eq. (B.7) 

A dummy index 
A dummy index 

See Eqs. (42), (43) and (44) 

See Eq. (62) 

An index for time step 
2 

Newton/meter 
Pressure 
Center pressure 

Hertzian pressure 

Radius of equivalent cylinder 
Rad iu s of cy 1 ind e r 1 

Radius of cylinder 2 

See Eq. (B.10) 

time 
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a 


a 

— a 



o 


0 



o 


Approach velocity 

center approach velocity 

Deformation velocity 
coordinate along film 

Coordinate separating the inlet and middle region 

Load per unit width of cylinder 
Dummy coordinate along film 
Pressure-viscosity coefficient 

Second pressure-viscosity coefficient 
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v-v 


M- 

u 



▼ (P) 

m 

A-Y (p) 
m r 


P 


P 


s 


P = 


P 


See Eq. (60) 

viscosity 
Ambient viscosity 

Poisson's ratio of cylinder 1 
Poisson's ratio of cylinder 2 
See Eqs. (56), (57) and (58) 

System equation 

Derivative of Y (p) with respect to p 

m i 

Density 

Ambient density 
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(b) 


Fig. 1-1 Geometry of the normal approach elas tohydrodynamic problem. 



Fig. 1-2 The relation between viscosity and pressure for a composite 
exponential lubricant. 
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14.0 x 10 


_ . Center Film 
1116 Thickness \ 

1 11.6 x 10 -5 \ 

2 8.0 x 10’ 5 
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\ Film Thickness 

\ 


5/ 4 


\\ ' 
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X = x/a 

Fig. 1-5 Pressure and deformation profiles, straight exponential 
lubricant, G = 3180, p = 10^ psi 



I 



lubricant, G = 3180, p = 1.25 x 10^ psi. 
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Fig. 1-12 Pressure and deformation profiles, composite exponential 
lubricant, p = 1.25 x 10^ psi, G = 3180. 



Fig. 1-13 Pressure and deformation profiles, composite exponential 
, p = 1.5 x 10 5 psi, G = 3180. 


lubricant 



















